以下のサイトからRとRStudioのインストーラをダウンロードし、インストール
R: https://cran.ism.ac.jp/ または https://cran.r-project.org/mirrors.html
RStudio: https://posit.co/download/rstudio-desktop/
左上のSource paneでR scriptにコードを記述する。
コードを選択してCmd+Return またはOption+Return
でコードが実行される。
Cmd+Option+P
で直前に実行したコードを再実行できる、コードを一部変更して確認する際に便利。
右上のEnvironment
paneには作成したオブジェクトが表示される。オブジェクトの中身を表示させたり、オブジェクトの削除も可能。
右下のFile
paneではフォルダ内のファイル、作成したグラフ、機能の説明などが表示される。Packagesパネルではパッケージのインストールとロードも可能。
基本的に1件のデータ解析ごとに専用のフォルダを作成し、それをR
projectファイルと関連させて管理する。
library(pacman)
p_load(tidyverse, readxl, openxlsx) # 一度に複数のパッケージをロードできる
Rは動的プログラミング言語(インタプリタ型言語)であり、Rコードを実行するためにコンパイルする必要は無く、1行ずつコードが解釈されて実行される。
四則演算は+, -, *, /で可能。ハッシュタグ
#はコードにおけるコメントを表す。
1 + 2 # #記号より右は無視される
## [1] 3
4 - 1
## [1] 3
12 * 6
## [1] 72
5 ^ 3 # べき乗
## [1] 125
20 / 3
## [1] 6.666667
20 %/% 3 # 整数商
## [1] 6
20 %% 3 # 剰余(mod)
## [1] 2
コロン:は1刻みで連続な数を簡単に作成できる。他の四則演算と組み合わせることも可能。
1:5
## [1] 1 2 3 4 5
-2:3
## [1] -2 -1 0 1 2 3
-1:-10
## [1] -1 -2 -3 -4 -5 -6 -7 -8 -9 -10
1:5 + 2
## [1] 3 4 5 6 7
1:5 - 2
## [1] -1 0 1 2 3
1:5 / 2
## [1] 0.5 1.0 1.5 2.0 2.5
1:5 * 2
## [1] 2 4 6 8 10
1:5 * 1:5
## [1] 1 4 9 16 25
Rには基本的な数学操作を行う関数や、データの統計的情報を知るための集計関数が最初から用意されている。
関数に渡すデータは引数(ひきすう)と呼ばれる。データそのもの、Rオブジェクト、他の関数の結果、いずれも引数になり得る。
下記のround関数におけるdigitsのように、関数で指定する対象の引数を仮引数(parameter)と呼び、仮引数に実際に
渡される値のことを実引数(argument)と呼ぶ。誤解が生じない場合は2つを区別せずに引数とする。
round(3.1415, digits = 2) # 数値を丸める、digitsは小数点以下の桁数
## [1] 3.14
factorial(4) # 階乗を計算(4*3*2*1)
## [1] 24
factorial(round(3.14, digits = 0)) # 関数の結果が引数になる場合、最も内側の関数から外側に向かって関数が実行される。
## [1] 6
x <- 1:50
length(1:20) # データの個数
## [1] 20
max(x) # 最大値
## [1] 50
which.max(x) # 最大値の位置(インデックス)
## [1] 50
min(x) # 最小値
## [1] 1
which.min(x) # 最小値の位置
## [1] 1
mean(x) # 平均値
## [1] 25.5
median(x) # 中央値
## [1] 25.5
sd(x) # 標準偏差
## [1] 14.57738
var(x) # 分散
## [1] 212.5
range(x) # 範囲
## [1] 1 50
quantile(x) # 分位点
## 0% 25% 50% 75% 100%
## 1.00 13.25 25.50 37.75 50.00
summary(x) # 要約統計量
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 13.25 25.50 25.50 37.75 50.00
関数の引数の指定は省略可能。またargs()関数で関数の引数名とデフォルト値を調べることができる。
round(3.1415, 2)
## [1] 3.14
args(round)
## function (x, digits = 0)
## NULL
# 数字の桁数を揃える
pi
## [1] 3.141593
format(pi, nsmall = 8) # 少数点以下8で揃える
## [1] "3.14159265"
format(pi, digits = 2) # 数字2桁以下を切り捨てる
## [1] "3.1"
# 0埋め(zero-padding)
formatC(5:12, width = 2, flag = "0")
## [1] "05" "06" "07" "08" "09" "10" "11" "12"
連続な数を生成する関数
seq(1, 10, length = 5) #1〜10を5分割で
## [1] 1.00 3.25 5.50 7.75 10.00
seq(1, 10, by = 2) #1から2ずつ増やして10まで
## [1] 1 3 5 7 9
seq(as.Date("2016-01-01"), as.Date("2016-01-05"), by = "day") # 日付も生成できる
## [1] "2016-01-01" "2016-01-02" "2016-01-03" "2016-01-04" "2016-01-05"
rep(1:3, times = 2) # (1, 2, 3)を2回
## [1] 1 2 3 1 2 3
rep(1:3, length = 5) #(1, 2, 3)を5個になるまで繰り返す
## [1] 1 2 3 1 2
rep(1:3, each = 3) #(1, 2, 3)の要素を3個ずつ
## [1] 1 1 1 2 2 2 3 3 3
replicate(10, 1+1) ## replicateはコマンドを複数回実行して結果をベクトルとして保存
## [1] 2 2 2 2 2 2 2 2 2 2
無作為標本抽出のための関数
sample(x = 1:6, size = 4) # ベクトルxからsize個の要素を取り出す
## [1] 2 4 3 5
sample(x = 1:6, size = 4, replace = TRUE) # replace = TRUEは取り出した要素を元に戻す(要素の重複が許される)
## [1] 3 6 4 6
x <- rnorm(100, mean = 0, sd = 1) # 平均0, 標準偏差1の正規分布からランダムに値を100個抽出
割合を指定して無作為抽出するにはsample_frac()を使う。ただし入力はデータフレーム。
sample_frac(data.frame(x = 1:20), size = 0.2) # 1〜20から20%の要素を取り出す
## x
## 1 9
## 2 19
## 3 15
## 4 7
Rではオブジェクトの中にデータを格納して保存することができる。オブジェクトを作るには名前を選び、
代入演算子 <-を使用する。
使用中のオブジェクトはEnvironmental
paneに表示される。同じ名前のオブジェクトは上書きされる。
a <- 1 # RStudioではOption + (-)キーが代入演算子のショートカット
3 -> b # 向きが逆でも大丈夫
a
## [1] 1
b
## [1] 3
ls() # すでに使用しているオブジェクト名がわかる
## [1] "a" "b" "x"
assign("C", 7, envir = globalenv()) # assign関数は代入演算子 -> と同じ、ただし環境を指定できる
代入操作をカッコで括ると代入とコンソールでの表示を同時に行う。
(c <- 1 + 3)
## [1] 4
Rのオブジェクトにはbasic type, mode, class,の3つの型がある:
typeof()を使用する。typeof()はvector
objectに対して使用するとmode型を返す。mode()を使用する。class()を使用する。Rでは基本的にデータをベクトル単位で扱う。データをc()で囲むとベクトルオブジェクトになる。ベクトルの要素単位で計算を実行する。
v1 <- c(8,7,2,5,1,3)
v2 <- c(1:6)
v1
## [1] 8 7 2 5 1 3
v2
## [1] 1 2 3 4 5 6
v1 - 2
## [1] 6 5 0 3 -1 1
v1 / 2
## [1] 4.0 3.5 1.0 2.5 0.5 1.5
v1 - v2
## [1] 7 5 -1 1 -4 -3
v1 * v2
## [1] 8 14 6 20 5 18
ベクトルは1つの同じ型のデータしか格納できない。
int <- c(1L, 5L) # 数字の後にLを付けると整数integerになる
text <- c("ace", "hearts")
int
## [1] 1 5
typeof(int)
## [1] "integer"
text
## [1] "ace" "hearts"
typeof(text)
## [1] "character"
行列は数学の行列とほぼ同じ、行と列からなる2次元の数値の配列。
die <- 1:6
m <- matrix(die, nrow = 2) ## 行列を作成、次元属性が付与される
m
## [,1] [,2] [,3]
## [1,] 1 3 5
## [2,] 2 4 6
m <- matrix(die, nrow = 2, byrow = T)
m
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
dimnames(m) <- list(c("r1","r2"), c("c1","c2","c3")) ## dimnames属性で行列の行と列に名前をつけることもできる
m
## c1 c2 c3
## r1 1 2 3
## r2 4 5 6
配列(array)は行列の拡張版であり、行列は配列の特殊な形。
配列は同じサイズの行列を数枚重ねたものであり、行と列以外に「層」を持つ。
ar <- array(c(11:14, 21:24, 31:34), dim = c(2,2,3)) ## array関数はn次元配列を作成, dim属性を指定する
ar
## , , 1
##
## [,1] [,2]
## [1,] 11 13
## [2,] 12 14
##
## , , 2
##
## [,1] [,2]
## [1,] 21 23
## [2,] 22 24
##
## , , 3
##
## [,1] [,2]
## [1,] 31 33
## [2,] 32 34
3行4列の行列を4層重ねて配列を作成。
Mat1 <- matrix(sample(1:12, 12, replace = TRUE), byrow = TRUE, nrow = 3)
Mat2 <- matrix(sample(1:12, 12, replace = TRUE), byrow = TRUE, nrow = 3)
Mat3 <- matrix(sample(1:12, 12, replace = TRUE), byrow = TRUE, nrow = 3)
Mat4 <- matrix(sample(1:12, 12, replace = TRUE), byrow = TRUE, nrow = 3)
Array1 <- array(c(Mat1, Mat2, Mat3, Mat4), dim = c(3, 4, 4))
class(Array1)
## [1] "array"
Array1
## , , 1
##
## [,1] [,2] [,3] [,4]
## [1,] 2 5 7 11
## [2,] 2 5 6 6
## [3,] 4 8 11 9
##
## , , 2
##
## [,1] [,2] [,3] [,4]
## [1,] 1 7 2 11
## [2,] 3 10 3 4
## [3,] 2 7 4 1
##
## , , 3
##
## [,1] [,2] [,3] [,4]
## [1,] 10 5 5 4
## [2,] 10 6 10 11
## [3,] 6 11 3 11
##
## , , 4
##
## [,1] [,2] [,3] [,4]
## [1,] 7 7 10 10
## [2,] 4 8 8 10
## [3,] 1 4 4 1
配列名[行番号, 列番号, 層番号]で要素を指定。
Array1[3, 1, 2] # 2番目の行列の3行1列目の要素を抽出
## [1] 2
Array1[ , , 3] # 3番目の行列を抽出
## [,1] [,2] [,3] [,4]
## [1,] 10 5 5 4
## [2,] 10 6 10 11
## [3,] 6 11 3 11
Array1[1:2, 1:2, ] # 各層から1-2行目と1-2列目の要素からなる2×2行列を抽出
## , , 1
##
## [,1] [,2]
## [1,] 2 5
## [2,] 2 5
##
## , , 2
##
## [,1] [,2]
## [1,] 1 7
## [2,] 3 10
##
## , , 3
##
## [,1] [,2]
## [1,] 10 5
## [2,] 10 6
##
## , , 4
##
## [,1] [,2]
## [1,] 7 7
## [2,] 4 8
Array1[2, , ] # 全ての行列の2行目を抽出、結果は行列になることに注意
## [,1] [,2] [,3] [,4]
## [1,] 2 3 10 4
## [2,] 5 10 6 8
## [3,] 6 3 10 8
## [4,] 6 4 11 10
リストは任意のRオブジェクトを1次元のグループにまとめる。
list1 <- list(100:130, "R", list(T, F))
names(list1) <- c("sequence", "character", "logical") ## namesで名前属性をつけることができる
list1
## $sequence
## [1] 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## [20] 119 120 121 122 123 124 125 126 127 128 129 130
##
## $character
## [1] "R"
##
## $logical
## $logical[[1]]
## [1] TRUE
##
## $logical[[2]]
## [1] FALSE
str(list1) # まとめているオブジェクトの型はstr()関数で見ることができる
## List of 3
## $ sequence : int [1:31] 100 101 102 103 104 105 106 107 108 109 ...
## $ character: chr "R"
## $ logical :List of 2
## ..$ : logi TRUE
## ..$ : logi FALSE
list1[1] # 番号で要素を選択、リスト形式で返される
## $sequence
## [1] 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## [20] 119 120 121 122 123 124 125 126 127 128 129 130
list_len3 <- vector("list", length = 3) # vector()で任意の長さのリストを作成できる
データフレームはリストの2次元バージョンであり、R版のExcelスプレッドシート。データフレームはベクトルをまとめて2次元の表にし、個々のベクトルは列になる。 それぞれのベクトルは異なる型のデータを格納できるが、同じ長さにする必要がある。言い換えるとデータフレームとは長さの等しいベクトルの名前付きリスト。データ分析のためには最も便利なストレージ構造。
df <- data.frame(face=c("ace","two","six"), suit=c("clubs","hearts","diamonds"), # names, class, row.names属性をもつ
value=c(1,2,3), stringsAsFactors = F)
df
## face suit value
## 1 ace clubs 1
## 2 two hearts 2
## 3 six diamonds 3
str(df)
## 'data.frame': 3 obs. of 3 variables:
## $ face : chr "ace" "two" "six"
## $ suit : chr "clubs" "hearts" "diamonds"
## $ value: num 1 2 3
typeof(df)
## [1] "list"
class(df)
## [1] "data.frame"
# dput()関数はオブジェクトをテキスト形式で出力する。defaultはコンソール上に出力される。
dput(df) # 出力結果をコピーしてRオブジェクトに代入すると、元と全く同一のオブジェクトが生成する。
## structure(list(face = c("ace", "two", "six"), suit = c("clubs",
## "hearts", "diamonds"), value = c(1, 2, 3)), class = "data.frame", row.names = c(NA,
## -3L))
tibbleはデータフレームの進化版。データフレームとは異なり、任意のRオブジェクト(データフレームやグラフを含む)を格納できる。 また入力データの型を自動変換しない。tidyverseではtibbleを標準の形式としている。 古い関数ではtibbleが使えないことがあるのでその場合はデータフレームを使用する。
as_tibble(iris) # as_tibble()でdataframeからtibbleへ変換できる
## # A tibble: 150 × 5
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## <dbl> <dbl> <dbl> <dbl> <fct>
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosa
## # ℹ 140 more rows
# as.data.frame()でさらにデータフレームに戻す、データフレームは全ての行を表示してしまうので、head()を使用
head(as.data.frame(as_tibble(iris)))
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
tibbleの例、tibbleはdefaultで上から10行までしか表示しない。
library(tidyverse)
library(nycflights13)
data("flights")
flights
## # A tibble: 336,776 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 -1 1004 1022
## 5 2013 1 1 554 600 -6 812 837
## 6 2013 1 1 554 558 -4 740 728
## 7 2013 1 1 555 600 -5 913 854
## 8 2013 1 1 557 600 -3 709 723
## 9 2013 1 1 557 600 -3 838 846
## 10 2013 1 1 558 600 -2 753 745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
# intは整数を表す
# dblはdoubleの略で実数を表す
# chrは文字ベクトルすなわち文字列を表す
# dttmはdate-times(日付+時刻)の略
# lglはlogicalの略で、TRUEまたはFALSEからなるベクトルを示す
# fctはfactorの略で、規定の値のカテゴリ変数をRで扱う時に使用
# dateは日付を表す
tibbleの列にリストが格納されている場合、その列を特にリスト列と呼ぶ。データフレーム形式ではリストは列として扱われ、含まれるベクトルの長さが異なるとエラーになる。tibble形式ではそのような制約は無い。
data.frame(x = list(1:3, 3:5))
## x.1.3 x.3.5
## 1 1 3
## 2 2 4
## 3 3 5
tibble(
x = list(1:3, 3:5, 5:10),
y = c("1, 2", "3, 4, 5", NA)
)
## # A tibble: 3 × 2
## x y
## <list> <chr>
## 1 <int [3]> 1, 2
## 2 <int [3]> 3, 4, 5
## 3 <int [6]> <NA>
角括弧[ ]内に添え字を記述して指定。添字は正の整数、負の整数、ゼロ、論理値、名前が使用できる。
vec_num <- c(1:30)
vec_num[7] # 10番目の要素、Rでは1から添え字が始まる
## [1] 7
vec_num[-3] # 3番目の要素以外
## [1] 1 2 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## [26] 27 28 29 30
vec_num[5:13] # :で範囲を指定
## [1] 5 6 7 8 9 10 11 12 13
vec_num[-(3:10)] # 3〜10番目以外
## [1] 1 2 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
vec_num[-3:-10] # こう書いても同じ
## [1] 1 2 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
vec_fruits <- c(1, 2, 3)
names(vec_fruits) <- c("apple", "banana", "orange") # 各要素の名前を付ける
vec_fruits["banana"] # 名前を指定して要素を抽出
## banana
## 2
ベクトルと同様、角括弧[ ]を使用。カンマで区切られた2つの添え字を記述して指定。
deck <- read.csv("deck.csv", header =T, sep = ",", stringsAsFactors = F) # csvファイルから読み込み、下記参照。
deck
## face suit value
## 1 king spades 13
## 2 queen spades 12
## 3 jack spades 11
## 4 ten spades 10
## 5 nine spades 9
## 6 eight spades 8
## 7 seven spades 7
## 8 six spades 6
## 9 five spades 5
## 10 four spades 4
## 11 three spades 3
## 12 two spades 2
## 13 ace spades 1
## 14 king clubs 13
## 15 queen clubs 12
## 16 jack clubs 11
## 17 ten clubs 10
## 18 nine clubs 9
## 19 eight clubs 8
## 20 seven clubs 7
## 21 six clubs 6
## 22 five clubs 5
## 23 four clubs 4
## 24 three clubs 3
## 25 two clubs 2
## 26 ace clubs 1
## 27 king diamonds 13
## 28 queen diamonds 12
## 29 jack diamonds 11
## 30 ten diamonds 10
## 31 nine diamonds 9
## 32 eight diamonds 8
## 33 seven diamonds 7
## 34 six diamonds 6
## 35 five diamonds 5
## 36 four diamonds 4
## 37 three diamonds 3
## 38 two diamonds 2
## 39 ace diamonds 1
## 40 king hearts 13
## 41 queen hearts 12
## 42 jack hearts 11
## 43 ten hearts 10
## 44 nine hearts 9
## 45 eight hearts 8
## 46 seven hearts 7
## 47 six hearts 6
## 48 five hearts 5
## 49 four hearts 4
## 50 three hearts 3
## 51 two hearts 2
## 52 ace hearts 1
deck[1, 2] # 1行2列目の要素を返す
## [1] "spades"
deck[1, 1:3] # 1行1〜3列、1:3をc(1,2,3)にしても同じ
## face suit value
## 1 king spades 13
deck[c(2, 6:11), 2:3] # 2行目と6〜11行目、2〜3列
## suit value
## 2 spades 12
## 6 spades 8
## 7 spades 7
## 8 spades 6
## 9 spades 5
## 10 spades 4
## 11 spades 3
deck[-(2:40), 1:3] # 2〜40行目以外
## face suit value
## 1 king spades 13
## 41 queen hearts 12
## 42 jack hearts 11
## 43 ten hearts 10
## 44 nine hearts 9
## 45 eight hearts 8
## 46 seven hearts 7
## 47 six hearts 6
## 48 five hearts 5
## 49 four hearts 4
## 50 three hearts 3
## 51 two hearts 2
## 52 ace hearts 1
deck[1, ] # スペースで列全体を取り出す
## face suit value
## 1 king spades 13
deck[ , 2] # スペースで行全体を取り出す
## [1] "spades" "spades" "spades" "spades" "spades" "spades"
## [7] "spades" "spades" "spades" "spades" "spades" "spades"
## [13] "spades" "clubs" "clubs" "clubs" "clubs" "clubs"
## [19] "clubs" "clubs" "clubs" "clubs" "clubs" "clubs"
## [25] "clubs" "clubs" "diamonds" "diamonds" "diamonds" "diamonds"
## [31] "diamonds" "diamonds" "diamonds" "diamonds" "diamonds" "diamonds"
## [37] "diamonds" "diamonds" "diamonds" "hearts" "hearts" "hearts"
## [43] "hearts" "hearts" "hearts" "hearts" "hearts" "hearts"
## [49] "hearts" "hearts" "hearts" "hearts"
deck[ , "suit"] # オブジェクトが名前を持っているならば名前で指定可能
## [1] "spades" "spades" "spades" "spades" "spades" "spades"
## [7] "spades" "spades" "spades" "spades" "spades" "spades"
## [13] "spades" "clubs" "clubs" "clubs" "clubs" "clubs"
## [19] "clubs" "clubs" "clubs" "clubs" "clubs" "clubs"
## [25] "clubs" "clubs" "diamonds" "diamonds" "diamonds" "diamonds"
## [31] "diamonds" "diamonds" "diamonds" "diamonds" "diamonds" "diamonds"
## [37] "diamonds" "diamonds" "diamonds" "hearts" "hearts" "hearts"
## [43] "hearts" "hearts" "hearts" "hearts" "hearts" "hearts"
## [49] "hearts" "hearts" "hearts" "hearts"
deck[16, "face"]
## [1] "jack"
deck[1:10, 1, drop=F] # 1列だけ抽出する場合、データフレームがほしい場合はdrop =Fとする
## face
## 1 king
## 2 queen
## 3 jack
## 4 ten
## 5 nine
## 6 eight
## 7 seven
## 8 six
## 9 five
## 10 four
class(deck[1:10, 1]) # 普通はベクトル型になる
## [1] "character"
ドル記号$はデータフレームやリストから名前の付いている要素をベクトルやデータフレームとして取り出すことができる。
deck$face # $でデータフレームから列を抽出
## [1] "king" "queen" "jack" "ten" "nine" "eight" "seven" "six" "five"
## [10] "four" "three" "two" "ace" "king" "queen" "jack" "ten" "nine"
## [19] "eight" "seven" "six" "five" "four" "three" "two" "ace" "king"
## [28] "queen" "jack" "ten" "nine" "eight" "seven" "six" "five" "four"
## [37] "three" "two" "ace" "king" "queen" "jack" "ten" "nine" "eight"
## [46] "seven" "six" "five" "four" "three" "two" "ace"
typeof(deck$value)
## [1] "integer"
list1 <- list(100:130, "R", list(T, F))
names(list1) <- c("sequence", "character", "logical") ## namesで名前属性をつけることができる
list1$sequence # ベクトル形式で取り出される
## [1] 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## [20] 119 120 121 122 123 124 125 126 127 128 129 130
typeof(list1$sequence)
## [1] "integer"
typeof(list1[1]) # 要素番号で取り出すとリスト形式になる
## [1] "list"
重要:リストの要素に名前が付いていない場合、2重の角括弧で要素を抽出すると、$記法と同じくベクトル/データフレーム形式になる
データフレーム形式になると角括弧[ ]を使用して、[行,
列]で要素が指定できる
list1[2]
## $character
## [1] "R"
typeof(list1[2])
## [1] "list"
list1[[2]]
## [1] "R"
typeof(list1[[2]])
## [1] "character"
list1$character
## [1] "R"
typeof(list1$character)
## [1] "character"
deck_list <- list(deck$face, deck$suit, deck)
deck_list[[3]][1,1,1] # データフレームなので[]で要素が指定できる
## [1] "king"
# deck_list[3][1,1,1] # deck_list[3]はリスト形式なのでエラーになる
attach()関数でデータフレームを読み込むと変数名だけで値にアクセスできるようになる。解除するにはdetach()を使う。
df <- data.frame(face=c("ace","two","six"), suit=c("spades","clubs","hearts"), # names, class, row.names属性をもつ
value=c(1,2,3), stringsAsFactors = F)
df$suit
## [1] "spades" "clubs" "hearts"
attach(df)
suit
## [1] "spades" "clubs" "hearts"
detach(df)
attach()と似た関数にwith()がある。構文はwith(data, expr, ...)。with関数内に限ってattach(data)した後と同様にexprが評価される。
df <- data.frame(face=c("queen","ace","six"), suit=c("diamonds","spades","clubs"), # names, class, row.names属性をもつ
value=c(1,2,3), stringsAsFactors = F)
with(df, paste(face, suit, sep = " of "))
## [1] "queen of diamonds" "ace of spades" "six of clubs"
式が真(TRUE)か偽(FALSE)かを判定するための演算子。式のそれぞれにTRUEかFALSEを返す。
ベクトルを比較する演算子を使った場合、Rは要素ごとに判定する。
1 < 2 # 左辺(LHS)は右辺(RHS)より小さい
## [1] TRUE
2 <= 2 # LHSはRHSと等しいか小さい
## [1] TRUE
5 > 6 # LHSはRHSより大きい
## [1] FALSE
4 >= 5 # LHSはRHSと等しいか大きい
## [1] FALSE
(2 + 3) == (4 + 1)
## [1] TRUE
((2 * 3) + 1) != (2 * (3 + 1))
## [1] TRUE
1 > 0:2
## [1] TRUE FALSE FALSE
c(1, 2, 3) == c(3, 2, 1)
## [1] FALSE TRUE FALSE
マッチ演算子 %in%
は左側の個々の値が右側のベクトルに含まれているかどうかをテストする。
従って左側の要素の個数だけTRUEまたはFALSEを返す。
1 %in% c(3, 4, 5)
## [1] FALSE
c(1, 2) %in% c(3, 4, 5)
## [1] FALSE FALSE
c(1, 2, 3, 4) %in% c(3, 4, 5)
## [1] FALSE FALSE TRUE TRUE
ブール演算子は複数の論理テストから1つのTRUEかFALSEを導出する。Rは個々の論理テストを実行してから、ブール演算子を使用して判定する。ベクトルとともに使うと、ブール演算子は要素単位で評価する。
a <- c(1, 2, 3)
b <- c(1, 2, 3)
c <- c(1, 2, 4)
a == b
## [1] TRUE TRUE TRUE
b == c
## [1] TRUE TRUE FALSE
!(b == c) # 論理テストの反転
## [1] FALSE FALSE TRUE
a == b & b == c # cond1 & cond2, cond1とcond2が両方ともTRUEかどうか
## [1] TRUE TRUE FALSE
a == b | b == c # cond1 & cond2, cond1とcond2のどちらかがTRUEかどうか
## [1] TRUE TRUE TRUE
any(a == b, b == c) # any(cond1, cond2, ...) condの中で1つでもTRUEになっているものがあるか
## [1] TRUE
all(a == b, b == c) # all(cond1, cond2,...) 全てのcondがTRUEかどうか
## [1] FALSE
ダブル演算子&&や||はそれぞれ&や|と同じように動作するが、ベクトル対応になっておらずスカラーに対して適用され、結果はTRUEまたはFALSE。
# a == b && b == c
# a == b || b == c
ベクトルの要素はTRUE/FALSEの配列で指定できる。
vec_num <- sample(1:6, 6, replace = FALSE) # 1:6をランダムに6個並べた配列を作る
vec_num
## [1] 2 1 3 6 4 5
vec_num[c(FALSE, TRUE, FALSE, TRUE, TRUE, FALSE)] # 2, 4, 5番目の要素を抽出
## [1] 1 6 4
従って論理テストの結果を利用して、要素を抽出することができる。sum()関数内でTRUEとFALSEはそれぞれ1と0に置換して計算されるので、条件に当てはまる要素の個数を数えることができる。また代入演算子を使用して値を置換することも可能。
vec_num2 <- sample(1:6, 20, replace = TRUE) # 1:6をランダムに20個並べた配列を作る
vec_num2
## [1] 6 3 3 1 2 3 4 3 2 2 6 5 5 4 3 6 5 6 6 2
(vec_num2 %% 2) == 0 # 要素が2で割り切れる(=偶数)かどうかを判定
## [1] TRUE FALSE FALSE FALSE TRUE FALSE TRUE FALSE TRUE TRUE TRUE FALSE
## [13] FALSE TRUE FALSE TRUE FALSE TRUE TRUE TRUE
vec_num2[(vec_num2 %% 2) == 0] # 偶数だけを抽出
## [1] 6 2 4 2 2 6 4 6 6 6 2
vec_num2 %in% c(2, 5) # 要素が2または5かどうかを判定
## [1] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE FALSE TRUE
## [13] TRUE FALSE FALSE FALSE TRUE FALSE FALSE TRUE
sum(vec_num2 %in% c(2, 5)) # 2または5の個数をカウント
## [1] 7
vec_num2[vec_num2 %in% c(2, 5)] # 2または5を抽出
## [1] 2 2 2 5 5 5 2
vec_num2[vec_num2 %in% c(2, 5)] <- 0 # 2または5を0に置換
vec_num2
## [1] 6 3 3 1 0 3 4 3 0 0 6 0 0 4 3 6 0 6 6 0
データフレームの場合も同様に論理テストが利用できる。
deck <- read.csv("deck.csv", header =T, sep = ",", stringsAsFactors = F) # csvファイルから読み込み
head(deck)
## face suit value
## 1 king spades 13
## 2 queen spades 12
## 3 jack spades 11
## 4 ten spades 10
## 5 nine spades 9
## 6 eight spades 8
sum(deck$face == "ace") # aceの個数を数える
## [1] 4
deck$value[deck$face == "queen"] # queenのポイントを表示
## [1] 12 12 12 12
deck[deck$suit == "spades", ] # spadesのカードを全て表示
## face suit value
## 1 king spades 13
## 2 queen spades 12
## 3 jack spades 11
## 4 ten spades 10
## 5 nine spades 9
## 6 eight spades 8
## 7 seven spades 7
## 8 six spades 6
## 9 five spades 5
## 10 four spades 4
## 11 three spades 3
## 12 two spades 2
## 13 ace spades 1
deck$suit == "hearts" & deck$face == "king" # ハートのキングを判定する論理テスト
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [37] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE FALSE FALSE
kingOfHearts <- deck$suit == "hearts" & deck$face == "king" # テストの結果をオブジェクトに格納
deck[kingOfHearts, ] # ハートのキングのカードを表示
## face suit value
## 40 king hearts 13
facecards <- deck$face %in% c("king", "queen", "jack") # フェイスカード(キング、クイーン、ジャック)を判定
deck[facecards, ]
## face suit value
## 1 king spades 13
## 2 queen spades 12
## 3 jack spades 11
## 14 king clubs 13
## 15 queen clubs 12
## 16 jack clubs 11
## 27 king diamonds 13
## 28 queen diamonds 12
## 29 jack diamonds 11
## 40 king hearts 13
## 41 queen hearts 12
## 42 jack hearts 11
deck$value[facecards]
## [1] 13 12 11 13 12 11 13 12 11 13 12 11
コードの塊を3回以上コピー&ペーストする場合、関数化を考えるべき。
関数を書くことのメリット
* コードが理解しやすい直感的な名前を関数につけることができる。 *
要求が変わっても、複数箇所を変更せずに、1ヶ所のコードを変更するだけで済む
* コピー&ペースト時の間違いをなくすことができる
df <- tibble::tibble(
a = rnorm(10),
b = rnorm(10),
c = rnorm(10),
d = rnorm(10)
)
df$a <- (df$a - min(df$a, na.rm = TRUE)) /
(max(df$a, na.rm = TRUE) - min(df$a, na.rm = TRUE))
df$b <- (df$b - min(df$b, na.rm = TRUE)) /
(max(df$b, na.rm = TRUE) - min(df$a, na.rm = TRUE)) ## 間違いあり
df$c <- (df$c - min(df$c, na.rm = TRUE)) /
(max(df$c, na.rm = TRUE) - min(df$c, na.rm = TRUE))
df$d <- (df$d - min(df$d, na.rm = TRUE)) /
(max(df$d, na.rm = TRUE) - min(df$d, na.rm = TRUE))
(df$a - min(df$a, na.rm = TRUE)) /
(max(df$a, na.rm = TRUE) - min(df$a, na.rm = TRUE))
## [1] 0.07320017 0.50851439 0.12038401 0.71248348 0.00000000 0.36848580
## [7] 1.00000000 0.64768651 0.79237052 0.18434257
x <- df$a
(x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
## [1] 0.07320017 0.50851439 0.12038401 0.71248348 0.00000000 0.36848580
## [7] 1.00000000 0.64768651 0.79237052 0.18434257
#> [1] 0.2892677 0.7509271 0.0000000 0.6781686 0.8530656 1.0000000 0.1716402
#> [8] 0.6107464 0.6116181 0.6008793
rng <- range(x, na.rm = TRUE) ## rangeは最小値と最大値を含んだベクトルを返す
(x - rng[1]) / (rng[2] - rng[1]) ## rng[1]は最小値、rng[2]は最大値
## [1] 0.07320017 0.50851439 0.12038401 0.71248348 0.00000000 0.36848580
## [7] 1.00000000 0.64768651 0.79237052 0.18434257
関数化のステップ:1. 関数に名前をつける, 2. 引数を関数の内側に書く, 3.
関数直後の{}にコードを書く。
rescale01 <- function(x) {
rng <- range(x, na.rm = TRUE)
(x - rng[1]) / (rng[2] - rng[1])
}
rescale01(c(0, 5, 10)) # 他の値で上手くいくか確認
## [1] 0.0 0.5 1.0
ベクトルの各要素に名前があるかどうかを示す論理ベクトルを返す関数。戻り値を指定しない場合、最後に計算した値を返す。
has_name <- function(x){
nms <- names(x)
if (is.null(nms)){
rep(FALSE, length(x))
} else{
!is.na(nms) & nms != ""
}
}
|| (or)や&&
(and)を使用して複数の論理式を結合できる。
||は最初にTRUEと評価された要素でTRUEを返し、残りの要素は評価しない。
&&ではFALSEと評価した最初の要素で、FALSEを返す。
if文では|や&を使用してはいけない、これらは複数値に対するベクトル演算を行う(なのでfilter()で使用する)。
論理ベクトルはany()やall()を使用して、単一値にまとめることができる。
等しいかどうかのテスト ==
もベクトル演算なので注意する。
a=c(1,1,1,1,1); b=c(1,1,0,0,1)
a == b
## [1] TRUE TRUE FALSE FALSE TRUE
identical(a, b) ## identical は非常に厳格で常に単一のTかFを返し、型強制はしない
## [1] FALSE
複数のif文をつなげることができる:
if (this) {
# do that
} else if (that) {
# do something else
} else {
}
switch関数を使用すると、位置や名前に基づいて選択したコードを評価できる。
op <- "1"
x <- 10
y <- 20
ans <- switch(op,
"1" = x + y,
"2" = x - y,
"3" = x * y,
"4" = x / y,
stop("Only can use 1, 2, 3, and 4"))
ans
## [1] 30
一つは計算するデータ、もう一つは計算の詳細を制御。例えばlog()ではデータがx、詳細が対数の底base。一般にデータ引数が最初に来る、詳細引数は最後に来て、通常はデフォルト値が存在。
# 正規近似を使って、平均値の信頼区間を計算
mean_ci <- function(x, conf = 0.95) { ## 0.95がdefaultの値
se <- sd(x) / sqrt(length(x))
alpha <- 1 - conf
mean(x) + se * qnorm(c(alpha / 2, 1 - alpha / 2))
}
x <- runif(100)
mean_ci(x)
## [1] 0.4164526 0.5251495
#> [1] 0.4976111 0.6099594
mean_ci(x, conf = 0.99)
## [1] 0.3993751 0.5422271
重み付き要約統計量を計算する関数
wt_mean <- function(x, w) {
sum(x * w) / sum(w)
}
wt_var <- function(x, w) {
mu <- wt_mean(x, w)
sum(w * (x - mu) ^ 2) / sum(w)
}
wt_sd <- function(x, w) {
sqrt(wt_var(x, w))
}
wt_mean(1:6, 1:5) ## xとwのベクトルの長さが異なるとエラーになる
## Warning in x * w:
## 長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません
## [1] 4.066667
重要な前提条件はチェックして、正しくないとstopと共にエラーを返すようにする
wt_mean <- function(x, w) {
if (length(x) != length(w)) {
stop("`x` and `w` must be the same length", call. = FALSE) ## xとwの長さが異なる場合、stopとともにエラーメッセージを返す
}
sum(w * x) / sum(w)
}
stopifnot()は各引数がTRUEかどうかを調べ、そうでないと汎用的なエラーメッセージを返す
wt_mean <- function(x, w, na.rm = FALSE) {
stopifnot(is.logical(na.rm), length(na.rm) == 1)
stopifnot(length(x) == length(w))
if (na.rm) {
miss <- is.na(x) | is.na(w)
x <- x[!miss]
w <- w[!miss]
}
sum(w * x) / sum(w)
}
# wt_mean(1:6, 6:1, na.rm = "foo")
Rの多くの関数では入力する引数の数に制限は無い、例えばsumやstr_cなど。
sum(1,2,3,4,5,6,7,8,9,10)
## [1] 55
stringr::str_c("a","b","c","d","e","f","g")
## [1] "abcdefg"
特別な引数「…」は入力引数がいくつでもマッチでき、他の関数に送り出すことができる。
commas <- function(...) stringr::str_c(..., collapse = ", ")
commas(letters[1:10])
## [1] "a, b, c, d, e, f, g, h, i, j"
rule <- function(..., pad = "-") {
title <- paste0(...)
width <- getOption("width") - nchar(title) - 5 ## getOption("width")はコンソール画面の幅
cat(title, " ", stringr::str_dup(pad, width), "\n", sep = "")
}
rule("Important output")
## Important output -----------------------------------------------------------
関数が返す値は通常は評価の最後の文。return()を使って早く返すように選択できる。考慮するべきは1. 値を早く返すと関数が読みやすくなるか 2. パイプを使って関数が書けるか。 明示的return文にするべきは主に次の3つのケース: 1. 入力が空の場合、2. 複雑なブロックが一つと簡単なブロックが一つある場合、3. パイプにできる関数を書く場合
重複をなくすツールの1つが関数化で、繰り返し使われるコードパターンを独立させて、たやすく再利用や更新できるようにする。もう1つのツールがイテレーション(iteration)で、複数の入力(異なる列やデーセット)に対して同じことを行う。イテレーションには命令型(for loop, while loop)と関数型がある。
forループの3つの要素:
出力:
ループを開始する前に、出力用のスペースを確保しておく必要がある。指定長の空ベクトルを作るにはvector()を使う。
シーケンス:
何でループするのかを決める。seq_along()は1からベクトルの長さ(行数)までの数列を生成、ベクトルの長さが0の時にもうまく処理してくれるので1:length()よりも安全。
本体: 実際の仕事をするコード。
df <- tibble(
a = rnorm(10),
b = rnorm(10),
c = rnorm(10),
d = rnorm(10)
)
# for loopで計算
output <- vector("double", ncol(df)) ## 1. 出力、結果を出力するための指定長の空ベクトルを作っておく
for (i in seq_along(df)){ ## 2. シーケンス
output[[i]] <- median(df[[i]]) ## 3. 本体
}
output
## [1] -0.0962729 0.3255362 -0.4405117 -0.6236106
[[]]は単一要素を扱うことが明示されるので、for loopの中では[[]]を使用する(例えアトミックベクトルでも[[]]の使用を推奨)。
forループには以下の4つのバリエーションがある:
出力は既にある、すなわち入力と同じ。データフレームを列のリストと考え、列をiteration。
df <- tibble(
a = rnorm(10),
b = rnorm(10),
c = rnorm(10),
d = rnorm(10)
)
# 関数で解く場合
rescale01 <- function(x) {
rng <- range(x, na.rm = TRUE)
(x - rng[1]) / (rng[2] - rng[1])
}
# for loopで解く場合
for (i in seq_along(df)){
df[[i]] <- rescale01(df[[i]]) ## 本体:rescale01()する
}
ベクトルでループするには基本的に3つの方法がある:
(i in seq_along(xs))で数値の添字を使って値xs[[i]]を抽出。names(xs)[[i]]で名前も取得できる。for (x in xs)、出力を効率的に格納するのは難しい。グラフを書いたりするなど、副作用だけを使う際に有効。(nm in names(xs))、xs[[nm]]で値にアクセスできる名前を指定。名前付きの出力を作るなら、結果ベクトルに名前を付けるのを忘れないようにする。results <- vector("list", length(x))
names(results) <- names(x)
出力の長さがどれだけになるかわからない場合もある。例:長さがランダムなベクトルをシミュレーションしたい。
means <- c(0, 1, 2)
output <- double()
for (i in seq_along(means)) {
n <- sample(100, 1) ## 1〜100の間の整数を一つランダムに選ぶ
#ベクトルoutputの最後にrnorm(..)で選んだn個の数字を追加して、ベクトルを長くしていく
output <- c(output, rnorm(n, means[[i]]))
}
str(output)
## num [1:152] -2.255 -0.864 0.621 0.143 0.951 ...
これでは毎回Rが前のデータを全てコピーしないといけないので、あまり効率的ではない。 より良い解法では、結果をリストに格納し、ループが終わってから一つのベクトルにまとめる。
out <- vector("list", length(means)) ## 長さが3のリストを作成できる
for (i in seq_along(means)) {
n <- sample(100, 1)
out[[i]] <- rnorm(n, means[[i]]) ## リストのi番目の要素に格納
}
str(out)
## List of 3
## $ : num [1:10] -0.594 1.064 1.116 -0.397 1.727 ...
## $ : num [1:93] 1.5048 1.8186 -0.0516 -0.3624 0.7088 ...
## $ : num [1:92] 3.31 2.88 3.91 2.02 1.27 ...
str(unlist(out)) ## ベクトルのリストを1つのベクトルにする、purrr::flatten_dbl()を使う方がベター
## num [1:195] -0.594 1.064 1.116 -0.397 1.727 ...
他の例:
paste(output, collapse="")でベクトルの要素を連結して文字列にする。dplyr::bind_rows(output)を使って連結して1つのデータフレームにまとめる。何回繰り返すか不明の場合はwhileを使用。例えば表が3回連続出るまでに何回硬貨投げを試行するのかを数えるwhileループ。
flip <- function() sample(c("T", "H"), 1)
flips <- 0
nheads <- 0
while (nheads < 3) {
if (flip() == "H") {
nheads <- nheads + 1
} else {
nheads <- 0
}
flips <- flips + 1
}
flips
## [1] 42
Rは関数型プログラミング言語であり、他の言語ほどforループが重要とは考えない。forループは関数でラップ(wrap)できる。
df <- tibble(
a = rnorm(10),
b = rnorm(10),
c = rnorm(10),
d = rnorm(10)
)
# for loopで各列の平均を計算
output <- vector("double", length(df))
for (i in seq_along(df)) {
output[[i]] <- mean(df[[i]])
}
output
## [1] -0.2607816 -0.8224295 0.1765180 -0.3967573
平均値の計算を関数化。中央値と標準偏差の計算も同様に関数化。
col_mean <- function(df) { ## 平均の計算を関数化
output <- vector("double", length(df))
for (i in seq_along(df)) {
output[i] <- mean(df[[i]])
}
output # 返り値
# return(output) # こう書いても同じ
}
col_mean(df)
## [1] -0.2607816 -0.8224295 0.1765180 -0.3967573
col_median <- function(df) { ## 中央値の計算を関数化
output <- vector("double", length(df))
for (i in seq_along(df)) {
output[i] <- median(df[[i]])
}
output
}
col_median(df)
## [1] -0.2797751 -0.5692534 0.2601980 -0.4112665
col_sd <- function(df) { #### 標準偏差の計算を関数化
output <- vector("double", length(df))
for (i in seq_along(df)) {
output[i] <- sd(df[[i]])
}
output
}
col_sd(df)
## [1] 0.4467616 1.0601581 0.4877669 1.0318189
平均値、中央値、標準偏差の関数を引数として渡せば、より一般的した関数を定義できる。関数を引数にする場合は()を付けない。
# 関数を別の関数に引き渡す
col_summary <- function(df, fun) { ## 引数funを追加、mean, median, sd等を入力すれば良い
out <- vector("double", length(df))
for (i in seq_along(df)) {
out[i] <- fun(df[[i]])
}
out
}
col_summary(df, median)
## [1] -0.2797751 -0.5692534 0.2601980 -0.4112665
col_summary(df, mean)
## [1] -0.2607816 -0.8224295 0.1765180 -0.3967573
purrrパッケージのマップ関数とは第1引数のベクトルの各要素について何かを行い、結果を保存するための関数。どの関数もベクトルを入力し、各要素に関数を適用し、入力と同じ長さ(同じ名前)の新たなベクトルを返す。第2引数.fは式、文字ベクトル、または整数ベクトルでも良い。3ドット「…」で.fが呼ばれる際の引数を渡す。文字や整数の場合はその位置の要素が抽出される(後述):
map(): makes a list., lapplyと同様map_lgl(): makes a logical vector.map_int(): makes an integer vector.map_dbl(): makes a double vector.map_chr(): makes a character vector.map(df, mean) # リストを返す、各列の名前も保持される
## $a
## [1] -0.2607816
##
## $b
## [1] -0.8224295
##
## $c
## [1] 0.176518
##
## $d
## [1] -0.3967573
map_dbl(df, mean) # ベクトルを返す
## a b c d
## -0.2607816 -0.8224295 0.1765180 -0.3967573
map_chr(df, mean) # 要素が文字型のベクトルを返す
## Warning: Automatic coercion from double to character was deprecated in purrr 1.0.0.
## ℹ Please use an explicit call to `as.character()` within `map_chr()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## a b c d
## "-0.260782" "-0.822429" "0.176518" "-0.396757"
map_dbl(df, median)
## a b c d
## -0.2797751 -0.5692534 0.2601980 -0.4112665
map_dbl(df, sd)
## a b c d
## 0.4467616 1.0601581 0.4877669 1.0318189
df %>% map_dbl(max) # パイプ演算子も使える
## a b c d
## 0.5653385 0.5763737 0.7829316 1.1426596
mtcarsデータセットをシリンダーの値で3つに分け、同じ線形モデルを適用したい。
models <- mtcars %>%
split(.$cyl) %>%
map(function(df) lm(mpg ~ wt, data = df)) ## 匿名関数を定義して使用
匿名関数の代わりに、one-sided fomulaオブジェクト(~)を使用して次のように書ける。ここでピリオドは現在のリスト要素を参照する代名詞のようなもの。
models <- mtcars %>%
split(.$cyl) %>%
map(~ lm(mpg ~ wt, data = .))
R^2のような要約統計量を取得するにはsummary()を実行してから、r.squaredという成分を抽出する。
models %>%
map(summary) %>% ## summaryを実行してから、リストに格納されているR2値成分を抽出する
map_dbl(~ .$r.squared)
## 4 6 8
## 0.5086326 0.4645102 0.4229655
# 文字でも抽出できる、こちらの方が一般的
models %>%
map(summary) %>%
map_dbl("r.squared")
## 4 6 8
## 0.5086326 0.4645102 0.4229655
整数を使ってその位置の要素を抽出することもできる。
x <- list(list(1, 2, 3), list(4, 5, 6), list(7, 8, 9))
x %>% map_dbl(2)
## [1] 2 5 8
lapply()は基本的にmap()と同じ。map()の方がpurrrの他の関数との一貫性があり、.fのショートカットを使う。sapply()はlapply()のwrapperで、出力を自動的に簡素化。どんな出力になるかわからないので自作関数に使うには問題がある。vapply()は型を定義する引数を追加するので、sapply()の安全版。入力が面倒になる。vapply(df, is.numeric, logical(1))はmap_lgl(df, is.numeric)と等価。map関数はベクトルしか作成できないのに対し、vapply()は行列も作成できる。map関数の練習問題:
# a. mtcarsの各列の平均を計算する
data("mtcars")
mtcars %>% map_dbl(mean)
## mpg cyl disp hp drat wt qsec
## 20.090625 6.187500 230.721875 146.687500 3.596563 3.217250 17.848750
## vs am gear carb
## 0.437500 0.406250 3.687500 2.812500
# b. nycflights13::flightsの各列の型を決定する
data("flights")
flights %>% map_chr(typeof)
## year month day dep_time sched_dep_time
## "integer" "integer" "integer" "integer" "integer"
## dep_delay arr_time sched_arr_time arr_delay carrier
## "double" "integer" "integer" "double" "character"
## flight tailnum origin dest air_time
## "integer" "character" "character" "character" "double"
## distance hour minute time_hour
## "double" "double" "double" "double"
# c. irisの各列の一意な値の個数を計算する
data("iris")
iris %>% map(unique) %>% map(length)
## $Sepal.Length
## [1] 35
##
## $Sepal.Width
## [1] 23
##
## $Petal.Length
## [1] 43
##
## $Petal.Width
## [1] 22
##
## $Species
## [1] 3
# 以下のようにより簡潔に書ける
iris %>% map_int(n_distinct)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 35 23 43 22 3
# d. mu = -10, 0, 10, 100のそれぞれについて10個の正規乱数を生成する
map(c(-10,0,10,100), rnorm, n = 10) ## rnorm()の引数は後ろにつける
## [[1]]
## [1] -9.693572 -9.105398 -9.972139 -11.043093 -9.622109 -9.986116
## [7] -10.782066 -9.902077 -11.154659 -8.897776
##
## [[2]]
## [1] -1.70493586 -0.11161431 -0.29043668 0.44228369 -1.74003872 -1.31723386
## [7] 0.03207093 0.64853885 -0.36374951 1.66519574
##
## [[3]]
## [1] 10.192985 9.696792 8.768456 9.734266 8.451277 9.981873 10.212536
## [8] 9.970477 11.647382 10.735504
##
## [[4]]
## [1] 100.45297 102.56478 99.23805 98.68471 100.63526 99.39195 98.66007
## [8] 100.41030 100.65494 99.92367
# 匿名関数を使用する場合
map(c(-10,0,10,100), function(mu) rnorm(n = 10, mean = mu))
## [[1]]
## [1] -10.692553 -11.468975 -9.678817 -8.837789 -10.170395 -10.416223
## [7] -8.806419 -9.002832 -9.165819 -11.525696
##
## [[2]]
## [1] -1.0843696 -0.4529417 -1.1100316 0.7267748 1.4531480 -1.2733976
## [7] -0.7647024 -0.2249176 1.8354420 0.4604920
##
## [[3]]
## [1] 9.633834 8.752625 9.148156 10.634914 11.378216 10.467256 9.944206
## [8] 8.456584 10.635121 10.217572
##
## [[4]]
## [1] 99.19464 100.77860 100.11078 97.18729 99.78555 99.02519 101.85769
## [8] 100.57560 99.61918 99.98137
# one-sided formulaオブジェクトを使用する場合
map(c(-10, 0, 10, 100), ~ rnorm(n = 10, mean = .))
## [[1]]
## [1] -9.182589 -9.737206 -9.824892 -8.893812 -9.167568 -9.988400
## [7] -10.590297 -9.541603 -10.193090 -12.055588
##
## [[2]]
## [1] -0.03976529 -0.40782951 1.04600315 0.72389716 -0.15182399 0.62638232
## [7] -0.58032425 1.00933751 -0.55416588 -0.41401237
##
## [[3]]
## [1] 10.349641 9.652933 10.646024 10.034130 10.263762 9.814422 11.122525
## [8] 9.144291 9.787309 8.294519
##
## [[4]]
## [1] 100.34735 99.94242 100.89763 99.04754 99.88458 101.51632 101.69118
## [8] 100.81106 100.69948 101.17728
# 3. map(1:5, runif)は何をするか、それはなぜか
map(1:5, runif)
## [[1]]
## [1] 0.6899524
##
## [[2]]
## [1] 0.1154215 0.5156386
##
## [[3]]
## [1] 0.9582620 0.5123572 0.5954985
##
## [[4]]
## [1] 0.31071324 0.99986653 0.72981389 0.02290406
##
## [[5]]
## [1] 0.1309487 0.1608921 0.6092839 0.3920713 0.1263547
# 以下と等価。
map(seq(1:5), runif)
## [[1]]
## [1] 0.1762757
##
## [[2]]
## [1] 0.3664654 0.7194777
##
## [[3]]
## [1] 0.8599465 0.1532712 0.5284390
##
## [[4]]
## [1] 0.20072910 0.18444203 0.06565291 0.73609323
##
## [[5]]
## [1] 0.7399966 0.4658728 0.7744339 0.3840219 0.5345716
# 4. map(-2:2, rnorm, n = 5)は何をするか、map_dbl(-2:2, rnorm, n = 5)はどうなるか
map(-2:2, rnorm, n = 5)
## [[1]]
## [1] -2.0593755 -1.6782743 -0.9160276 -2.3464813 0.2129120
##
## [[2]]
## [1] -0.5782569 -0.3170911 -2.6440907 -1.6742033 -1.3217819
##
## [[3]]
## [1] 0.85090377 0.65089650 0.09451189 0.10839872 1.13291065
##
## [[4]]
## [1] 1.211274 -0.739638 2.442857 2.185321 2.414876
##
## [[5]]
## [1] 2.488894 1.520707 3.071209 2.076521 1.454397
# map_dbl(-2:2, rnorm, n = 5) # エラーになる
関数safely()は関数を引数に取り、修正版を返す。この場合、修正された関数はエラーを投げない。代わりに次の2つのオブジェクトをリストにまとめて返す:
result: 元の結果。エラーがあればNULL。
error: エラーオブジェクト。処理が成功すればNULL。
関数が成功すれば、要素resultには結果が入り、要素errorはNULLになる。
safe_log <- safely(log) ## safelyはresultとerrorを1つのリストにまとめて返す
str(safe_log(10))
## List of 2
## $ result: num 2.3
## $ error : NULL
str(safe_log("a"))
## List of 2
## $ result: NULL
## $ error :List of 2
## ..$ message: chr " 数学関数に数値でない引数が渡されました "
## ..$ call : language .Primitive("log")(x, base)
## ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"
safely()はmapで使うことができる。list(result,
error)を要素とするリストが得られる。
x <- list(1, 10, "a")
y <- x %>% map(safely(log))
str(y)
## List of 3
## $ :List of 2
## ..$ result: num 0
## ..$ error : NULL
## $ :List of 2
## ..$ result: num 2.3
## ..$ error : NULL
## $ :List of 2
## ..$ result: NULL
## ..$ error :List of 2
## .. ..$ message: chr " 数学関数に数値でない引数が渡されました "
## .. ..$ call : language .Primitive("log")(x, base)
## .. ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"
purrr::list_transpose()を得られたリストに対して使うと、resultとerror、それぞれのリスト(list(result), list(error))に変換される。list_transpose()は行列転置関数t()のlist版。 例えば、pair of lists <=> list of pairs。 data.frameをlistとして渡すとlist of rowsが得られる。
y <- y %>% list_transpose()
str(y)
## List of 2
## $ result:List of 3
## ..$ : num 0
## ..$ : num 2.3
## ..$ : NULL
## $ error :List of 3
## ..$ : NULL
## ..$ : NULL
## ..$ :List of 2
## .. ..$ message: chr " 数学関数に数値でない引数が渡されました "
## .. ..$ call : language .Primitive("log")(x, base)
## .. ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"
possibly()はsafely()同様、常に成功する、エラーの際default値を返す。quietly()はエラーを捕捉せずに出力、メッセージ、警告を補足する。
x <- list(1, 10, "a")
x %>% map_dbl(possibly(log, NA_real_))
## [1] 0.000000 2.302585 NA
x <- list(1, -1)
x %>% map(quietly(log)) %>% str()
## List of 2
## $ :List of 4
## ..$ result : num 0
## ..$ output : chr ""
## ..$ warnings: chr(0)
## ..$ messages: chr(0)
## $ :List of 4
## ..$ result : num NaN
## ..$ output : chr ""
## ..$ warnings: chr " 計算結果が NaN になりました "
## ..$ messages: chr(0)
単一の入力ではなく、複数の関連する入力について並列にイテレーションする場合、関数map2()とpmap()を使う。以下は平均値の異なる正規乱数を複数シミュレーションする例。
# 平均値の異なる正規乱数を複数シミュレーションする
mu <- list(5, 10, -3) # 試す平均値のリスト
# mu <- c(5, 10, -3) ## ベクトルにしても結果は同じ
mu %>%
map(rnorm, n = 5) %>%
str()
## List of 3
## $ : num [1:5] 3.8 6.25 4.66 5.32 2.8
## $ : num [1:5] 9.38 10.44 10.27 9.18 11.62
## $ : num [1:5] -3.13 -4.09 -1.75 -3.11 -1.48
平均値に加えて標準偏差sigmaも変化させる場合。平均と標準偏差のベクトルに添え字を当ててiterationするやり方。
sigma <- list(1, 5, 10)
# sigma <- c(1, 5, 10)
seq_along(mu) %>%
map(~ rnorm(n = 5, mean = mu[[.]], sd = sigma[[.]])) %>% # n =, mean = , sd =, は省略可
str()
## List of 3
## $ : num [1:5] 5.37 2.84 5.35 5.93 5.74
## $ : num [1:5] 18.32 12.18 12.21 9.17 7.93
## $ : num [1:5] -4.833 -0.727 -5.323 0.499 1.683
引数を2つ取れる、map2()で2つのベクトルを並列に処理。コードの意図がわかりやすくなる。変化する引数は関数の前に、変化しない引数は関数の後に書くことに注意。
map2(mu, sigma, rnorm, n = 5) %>%
str()
## List of 3
## $ : num [1:5] 6.33 5.4 3.72 6.22 6.15
## $ : num [1:5] 13.46 14.33 8.86 11.68 8.82
## $ : num [1:5] -1.06 -4.58 6.67 -9.31 6.25
map3(), map4(),..等も考えられるが、purrrには代わりに引数のリストを引数として取るpmap()が用意されている。以下は平均値、標準偏差、サンプルサイズを変化させる場合。
mu <- list(5, 10, -3) # 試す平均値のリスト
sigma <- list(1, 5, 10)
n <- list(1, 3, 5)
args1 <- list(n, mu, sigma)
args1 %>%
pmap(rnorm) %>%
str()
## List of 3
## $ : num 4.48
## $ : num [1:3] 7.41 10.66 7.51
## $ : num [1:5] -19.399 7.398 4.153 0.328 -12.576
リストの要素に名前を付けないと、pmap()は位置を代わりに使うので、引数に名前を付けた方が良い。
args2 <- list(mean = mu, sd = sigma, n = n) ## 引数に名前を付けた方がわかりやすい
args2 %>%
pmap(rnorm) %>%
str()
## List of 3
## $ : num 4.72
## $ : num [1:3] 8.231 4.867 -0.321
## $ : num [1:5] -2.91 -11.074 -11.724 -10.815 -0.879
引数は全て同じ長さなので、データフレームに格納しておくとより便利でわかりやすい。
params <- tribble(
~mean, ~sd, ~n,
5, 1, 1,
10, 5, 3,
-3, 10, 5
)
params %>%
pmap(rnorm)
## [[1]]
## [1] 3.769903
##
## [[2]]
## [1] -0.7947732 9.6287552 6.5797799
##
## [[3]]
## [1] 10.3802031 -4.3383396 -0.4663692 -10.7885070 8.9599208
もう一段階複雑にして、関数への引数だけでなく、関数も変化させるにはinvoke_map()を使用する。第1引数は関数のリスト、第2引数は関数ごとに異なる引数を指定するリストのリスト、第3引数以降は全関数に渡される。
f <- c("runif", "rnorm", "rpois")
param <- list(
list(min = -1, max = 1),
list(sd = 5),
list(lambda = 10)
)
invoke_map(f, param, n = 5) %>% str()
## Warning: `invoke_map()` was deprecated in purrr 1.0.0.
## ℹ Please use map() + exec() instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## List of 3
## $ : num [1:5] 0.964 -0.384 0.238 0.489 0.565
## $ : num [1:5] -9.06 2.68 -7.17 -4.65 3.49
## $ : int [1:5] 8 9 9 10 15
パラメータと計算結果を1つのデータフレームにまとめることもできる。
sim <- tribble( ## データフレームに関数と引数のリストを格納
~f, ~params,
"runif", list(min = -1, max = 1),
"rnorm", list(sd = 5),
"rpois", list(lambda = 10)
)
sim %>%
mutate(sim = invoke_map(f, params, n = 10)) ## シミュレーション結果を列simに格納
## # A tibble: 3 × 3
## f params sim
## <chr> <list> <list>
## 1 runif <named list [2]> <dbl [10]>
## 2 rnorm <named list [1]> <dbl [10]>
## 3 rpois <named list [1]> <int [10]>
ウォークは関数をその戻り値ではなく、副作用のために呼び出したい際にmapの代わりに使う。関数を適用しつつ第一引数を暗黙的に返す。グラフのリストとファイル名のベクトルがあると、pwalk()を使って、各ファイルをディスクの対応する場所に保存できる。
plots <- mtcars %>%
split(.$cyl) %>%
map(~ ggplot(., aes(mpg, wt)) + geom_point())
paths <- stringr::str_c(names(plots), ".pdf")
pwalk(list(paths, plots), ggsave, path = tempdir())
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
purrrにはforループの他のパターンを抽象化する関数も用意されている。
TRUEかFALSEを返す述語関数とともに機能する関数、keep()とdiscard()は述語がTRUEまたはFALSEとなる要素をそれぞれ返す。
iris %>%
keep(is.factor) %>%
str()
## 'data.frame': 150 obs. of 1 variable:
## $ Species: Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
iris %>%
discard(is.factor) %>%
str()
## 'data.frame': 150 obs. of 4 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
some()とevery()は述語が要素のいずれか、またはすべてでTRUEとなるかどうかを判定する。
x <- list(1:5, letters, list(10))
x %>%
some(is_character)
## [1] TRUE
x %>%
every(is_vector)
## [1] TRUE
detect()は述語が真となる最初の要素を返す、detect_index()はその位置を返す。
x <- sample(10)
x
## [1] 9 3 7 8 5 4 10 6 2 1
x %>%
detect(~ . > 5)
## [1] 9
x %>%
detect_index(~ . > 5)
## [1] 1
head_while()とtail_while()は述語がTRUEである先頭または末尾からの要素列を返す。
x %>%
head_while(~ . > 5)
## [1] 9
x %>%
tail_while(~ . > 5)
## integer(0)
reduce()は「2項」関数(2つの主入力がある)をとり、単一要素が残るまでリストに繰り返しその関数を適用する。例えばデータフレームのリストに対して、要素を結合して単一データフレームにする場合に使用。
dfs <- list(
age = tibble(name = "John", age = 30),
sex = tibble(name = c("John", "Mary"), sex = c("M", "F")),
trt = tibble(name = "Mary", treatment = "A")
)
dfs %>% reduce(full_join)
## Joining with `by = join_by(name)`
## Joining with `by = join_by(name)`
## # A tibble: 2 × 4
## name age sex treatment
## <chr> <dbl> <chr> <chr>
## 1 John 30 M <NA>
## 2 Mary NA F A
あるいはベクトルのリストがあり、共通部分を求めたい場合。
vs <- list(
c(1, 3, 5, 6, 10),
c(1, 2, 3, 7, 8, 10),
c(1, 2, 3, 4, 8, 9, 10)
)
vs %>% reduce(intersect) ## intersect()は両方に含まれている要素を抽出する
## [1] 1 3 10
accumulate()もreduce()と同様に機能するが、中間結果を全て保存する。累積和を求めることができる。
x <- sample(10)
x
## [1] 9 8 1 6 4 2 10 3 7 5
x %>% accumulate(`+`) ## accumulateは中間結果を全て保存する
## [1] 9 17 18 24 28 30 40 43 50 55
Rにおける文字列の操作と正規表現(regular
expression)。stringrパッケージを使用する。
library(tidyverse)
library(stringr)
data("words")
特殊文字はバックスラッシュ\を使ってエスケープする。文字列の中身そのものを表示させるにはwriteLines()。
double_quote <- "\"" ## バックスラッシュでエスケープ
double_quote
## [1] "\""
writeLines(double_quote) ## 文字列の中身そのものを調べるときはwriteLines
## "
return <- "a \nb " ## 改行
tab <- "a \t b" ## タブ
mu <- "\u00b5"
writeLines(return)
## a
## b
writeLines(tab)
## a b
writeLines(mu)
## µ
?"'" ## 特殊文字一覧
文字列を操作するstringrの関数は全て「str_」で始まるので、RStudioの自動補完機能で一覧が表示される。
str_length(c("a", "R for data science", NA)) ## 文字列の長さ
## [1] 1 18 NA
str_c("x", "y") ## 文字列の連結
## [1] "xy"
str_c("x", "y", sep = ", " ) ## sep引数を使って区切り文字を制御
## [1] "x, y"
x <- c("abc", NA)
str_c("|-", x, "-|") ## NAはそのまま
## [1] "|-abc-|" NA
## "NA"と出力するにはstr_replace_naを使用
str_c("|-", str_replace_na(x), "-|")
## [1] "|-abc-|" "|-NA-|"
str_c()はベクトル化関数で、要素の少ないベクトルを最長ベクトルの要素数に自動的にリサイクルして合わせる。
str_c("prefix-", c("a", "b", "c"), "-suffix")
## [1] "prefix-a-suffix" "prefix-b-suffix" "prefix-c-suffix"
引数collapse = ""とすると、文字列のベクトルをまとめて一つの文字列にする。
str_c(c("x", "y", "z"), "a", collapse = "")
## [1] "xayaza"
str_c(c("x", "y", "z"), "a")
## [1] "xa" "ya" "za"
str_c(c("x", "y", "z"), collapse = "")
## [1] "xyz"
str_sub()は文字列が短すぎてもエラーにならず、できるだけ部分文字列を取り出す。
x <- c("Apple", "Banana", "Pear")
str_sub(x, 1, 3) ## 1〜3文字目を取り出す
## [1] "App" "Ban" "Pea"
str_sub(x, -3, -1) ## 負数は末尾からの位置
## [1] "ple" "ana" "ear"
# str_sub()の代入形式を使って、文字列を変更することもできる
str_sub(x ,1, 1) <- str_to_lower(str_sub(x, 1, 1))
x
## [1] "apple" "banana" "pear"
後述の正規表現とstr_extract()を使用することもできる。
x <- c("Apple", "Banana", "Pear")
str_extract(x, "^.{2}") # 先頭の2文字を取り出す
## [1] "Ap" "Ba" "Pe"
str_extract(x, ".{3}$") # 末尾の3文字を取り出す
## [1] "ple" "ana" "ear"
str_view()は文字列と正規表現を取り、マッチがどうなるかを示す。最も単純なパターンマッチは厳密なマッチ。改行以外のどんな文字ともマッチする「.」。
x <- c("apple", "banana", "pear")
str_view(x, "an") ## 厳密なマッチ
## [2] │ b<an><an>a
str_view(x, ".a.") ## .は任意の一文字
## [2] │ <ban>ana
## [3] │ p<ear>
正規表現の「.」ではなく、文字ピリオド「.」をマッチさせるには、エスケープして「\.」とする。ただしR上では文字列のエスケープと認識されてしまい、「\.」のような文字列エスケープは無いためエラーとなる。この問題を回避するため正規表現の\.を作るには、R上ではバックスラッシュを2つ重ねて「\\.」とする。参考:https://opur.club/textbook/2018-5-2/, https://www.jaysong.net/RBook/string.html#%E6%AD%A3%E8%A6%8F%E8%A1%A8%E7%8F%BE
# To create the regular expression, we need \\
dot <- "\\."
# But the expression itself only contains one:
writeLines(dot)
## \.
#> \.
# And this tells R to look for an explicit .
str_view(c("abc", "a.c", "bef"), "a\\.c")
## [2] │ <a.c>
文字通りのバックスラッシュとマッチさせるには、「\\\\」と、4つ重ねる必要がある。
x <- "a\\b"
writeLines(x)
## a\b
#> a\b
str_view(x, "\\\\")
## [1] │ a<\>b
| マッチ | 正規表現 | Rでの表現 |
|---|---|---|
| 任意の数字 | \d |
\\d |
| 空白文字(空白, タブ, 改行) | \s |
\\s |
| 記号を除く文字 | \w |
\\w |
| 単語の境界 | \b |
\\b |
ピリオド「.」 |
\. |
\\. |
バックスラッシュ「\」 |
\\ |
\\\\ |
x <- c("apple", "banana", "pear")
str_view(x, "^a")
## [1] │ <a>pple
str_view(x, "a$")
## [2] │ banan<a>
x <- c("apple pie", "apple", "apple cake")
str_view(x, "apple")
## [1] │ <apple> pie
## [2] │ <apple>
## [3] │ <apple> cake
## 文字列全体とだけマッチするには^と$の両方でアンカーする
str_view(x, "^apple$")
## [2] │ <apple>
| マッチ | 正規表現 |
|---|---|
| a, b, cのどれか1文字 | [abc] |
| a,b,c以外ならどの文字とも | [^abc] |
| この並びで現れる | (abc) |
# Look for a literal character that normally has special meaning in a regex
str_view(c("abc", "a.c", "a*c", "a c"), "a[.]c")
## [2] │ <a.c>
str_view(c("abc", "a.c", "a*c", "a c"), ".[*]c")
## [3] │ <a*c>
str_view(c("abc", "a.c", "a*c", "a c"), "a[ ]")
## [4] │ <a >c
# 複数の候補パターンから選ぶには候補指定|を使用する
str_view(c("grey", "gray"), "gr(e|a)y")
## [1] │ <grey>
## [2] │ <gray>
| 意味 | 正規表現 |
|---|---|
| 0か1 | ? |
| 1以上 | + |
| 0以上 | * |
| n個 | {n} |
| n個以上 | {n,} |
| n個以下 | {,n} |
| m個以上、n個以下 | {m,n} |
x <- "1888 is the longest year in Roman numerals: MDCCCLXXXVIII"
str_view(x, "CC?")
## [1] │ 1888 is the longest year in Roman numerals: MD<CC><C>LXXXVIII
str_view(x, "CC+")
## [1] │ 1888 is the longest year in Roman numerals: MD<CCC>LXXXVIII
str_view(x, 'C[LX]+') ## [LX]はLXのどれかとマッチする
## [1] │ 1888 is the longest year in Roman numerals: MDCC<CLXXX>VIII
str_view(x, "C{2}") ## 2回繰り返す
## [1] │ 1888 is the longest year in Roman numerals: MD<CC>CLXXXVIII
str_view(x, "C{2,}") ## 2回以上繰り返す
## [1] │ 1888 is the longest year in Roman numerals: MD<CCC>LXXXVIII
str_view(x, "C{2,3}") ## 2〜3回繰り返す
## [1] │ 1888 is the longest year in Roman numerals: MD<CCC>LXXXVIII
str_view(x, 'C{2,3}?') ## defaultでは最長文字列と最長マッチする、?を後ろにつけて最短文字列とマッチさせることもできる
## [1] │ 1888 is the longest year in Roman numerals: MD<CC>CLXXXVIII
str_view(x, 'C[LX]+?')
## [1] │ 1888 is the longest year in Roman numerals: MDCC<CL>XXXVIII
| 意味 | 正規表現 |
|---|---|
| 先頭がA〜Fで始まる | ^[A-F] |
| 末尾が数字の1〜9で終わる | [1-9]$ |
| 何らかの文字が0回以上出現 | .* |
| 3桁の数字 | [0-9]{3} |
| 3桁-3桁-4桁の電話番号 | ^[0-9]{3}-[0-9]{3}-[0-9]{4}$ |
()でグループ化すると、キャプチャも同時に行われる。キャプチャとは、()内の文字列を記憶し、後で参照できるようにする機能。記憶した文字列は、後で「\1」,
「\2」,
「\3」,…のように指定することで利用することができる。正規表現内の()の数だけグループができる。以下のwordsは980個の英単語のリスト。
## ()は番号付きのcapturing groupを作成、\1, \2によりそれらを後方参照できる
str_view(fruit, "(..)\\1", match = TRUE)
## [4] │ b<anan>a
## [20] │ <coco>nut
## [22] │ <cucu>mber
## [41] │ <juju>be
## [56] │ <papa>ya
## [73] │ s<alal> berry
str_view(fruit, "(.)(.)\\2\\1", match = TRUE)
## [5] │ bell p<eppe>r
## [17] │ chili p<eppe>r
str_view(words, "(.).\\1.\\1", match = TRUE)
## [265] │ <eleve>n
文字ベクトルがパターンにマッチするかどうかの確認は、str_detect()を使う。
x <- c("apple", "banana", "pear")
str_detect(x, "e") ## 入力と同じ長さの論理ベクトルを返す
## [1] TRUE FALSE TRUE
数値表現ではTRUEは1,
FALSEは0としてカウントされる。sum()とmean()でそれぞれ、個数と割合がわかる。
str_detect()とgrepl()はほぼ同じ、文法が異なる。
x <- c("apple", "banana", "pear")
grepl("e", x)
## [1] TRUE FALSE TRUE
# 一般単語でtで始まる単語の個数
sum(str_detect(words, "^t"))
## [1] 65
# 一般単語で母音で終わる単語の割合
# [aeiou]はa, e, i, o, uのいずれか、$は末尾を指定するアンカー
mean(str_detect(words, "[aeiou]$"))
## [1] 0.2765306
sum(str_detect(words, "[aeiou]$"))
## [1] 271
母音を含まない単語を探し出す2つの方法。前者の方がわかりやすい。
# 少なくとも母音を1つ含む単語をすべて探し出して、補集合を取る
no_vowels_1 <- !str_detect(words, "[aeiou]")
# 子音だけからなる単語を全て探し出す
# [^aeiou]はaeiou以外、^と$で囲むと文字列全体、+は一回以上の繰り返し
no_vowels_2 <- str_detect(words, "^[^aeiou]+$")
identical(no_vowels_1, no_vowels_2) ## identicalで比較
## [1] TRUE
str_detect()を使用してから、論理ベクトルでリストの部分集合を取る。またはstr_subse()を使う。
# 論理ベクトルで部分集合を取る
words[str_detect(words, "x$")]
## [1] "box" "sex" "six" "tax"
# str_subsetを使用
str_subset(words, "x$")
## [1] "box" "sex" "six" "tax"
通常は対象はデータフレームの列にあるので、filter()を使う。
df <- tibble(
word = words,
# seq_along(x)はベクトルxの長さと同じ数の数列を作成、1:length(x)と等価、インデックスを簡単に作れる
i = seq_along(word)
)
df %>%
filter(str_detect(word, "x$"))
## # A tibble: 4 × 2
## word i
## <chr> <int>
## 1 box 108
## 2 sex 747
## 3 six 772
## 4 tax 841
str_count()は単なるマッチの可否ではなく、文字列にマッチがいくつあるかを示す。
x <- c("apple", "banana", "pear")
str_count(x, "a") ## 文字列にマッチがいくつあるか
## [1] 1 3 1
# 平均すると、1つの単語にどれだけ母音があるか?
mean(str_count(words, "[aeiou]"))
## [1] 1.991837
df %>%
mutate(
vowels = str_count(word, "[aeiou]"),
consonants = str_count(word, "[^aeiou]")
)
## # A tibble: 980 × 4
## word i vowels consonants
## <chr> <int> <int> <int>
## 1 a 1 1 0
## 2 able 2 2 2
## 3 about 3 3 2
## 4 absolute 4 4 4
## 5 accept 5 2 4
## 6 account 6 3 4
## 7 achieve 7 4 3
## 8 across 8 2 4
## 9 act 9 1 2
## 10 active 10 3 3
## # ℹ 970 more rows
マッチが重なることは絶対に無いことに注意。
str_count("abababa", "aba")
## [1] 2
str_view_all("abababa", "aba")
## Warning: `str_view_all()` was deprecated in stringr 1.5.0.
## ℹ Please use `str_view()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [1] │ <aba>b<aba>
実際にマッチしたテキストを抽出するにはstr_extract()を使う。解析サンプルとしてstringr::sentencesに格納されているHarvard英語コーパスを使用。720個の文がリストに格納されている。
length(sentences)
## [1] 720
head(sentences)
## [1] "The birch canoe slid on the smooth planks."
## [2] "Glue the sheet to the dark blue background."
## [3] "It's easy to tell the depth of a well."
## [4] "These days a chicken leg is a rare dish."
## [5] "Rice is often served in round bowls."
## [6] "The juice of lemons makes fine punch."
色名を含む文を全て探す。最初に色名ベクトルを作成し、それを一つの正規表現にする。
colors <- c("red", "orange", "yellow", "green", "blue", "purple")
color_match <- str_c(colors, collapse = "|")
color_match
## [1] "red|orange|yellow|green|blue|purple"
色を含む文の部分集合リストを作り、どの色か調べるために色を抽出する。
has_color <- str_subset(sentences, color_match) ## 色を含む文を選択してリストを作る
matches <- str_extract(has_color, color_match) ## 色を抽出
# head(matches)
str_view_all(has_color, color_match)
## [1] │ Glue the sheet to the dark <blue> background.
## [2] │ Two <blue> fish swam in the tank.
## [3] │ The colt rea<red> and threw the tall rider.
## [4] │ The wide road shimme<red> in the hot sun.
## [5] │ See the cat glaring at the sca<red> mouse.
## [6] │ A wisp of cloud hung in the <blue> air.
## [7] │ Leaves turn brown and <yellow> in the fall.
## [8] │ He orde<red> peach pie with ice cream.
## [9] │ Pure b<red> poodles have curls.
## [10] │ The spot on the blotter was made by <green> ink.
## [11] │ Mud was spatte<red> on the front of his white shirt.
## [12] │ The sofa cushion is <red> and of light weight.
## [13] │ The sky that morning was clear and bright <blue>.
## [14] │ Torn scraps litte<red> the stone floor.
## [15] │ The doctor cu<red> him with these pills.
## [16] │ The new girl was fi<red> today at noon.
## [17] │ The third act was dull and ti<red> the players.
## [18] │ A <blue> crane is a tall wading bird.
## [19] │ Live wires should be kept cove<red>.
## [20] │ It is hard to erase <blue> or <red> ink.
## ... and 37 more
str_extract()は最初のマッチしか抽出しないことに注意。
# マッチを2つ以上含む文を選択
more <- sentences[str_count(sentences, color_match) >1]
str_view_all(more, color_match)
## [1] │ It is hard to erase <blue> or <red> ink.
## [2] │ The <green> light in the brown box flicke<red>.
## [3] │ The sky in the west is tinged with <orange> <red>.
str_extract(more, color_match)
## [1] "blue" "green" "orange"
マッチを全て取得するにはstr_extract_all()を使用する。返り値が2つ以上になるので、ベクトルではなくリストが返ってくる。
str_extract_all(more, color_match)
## [[1]]
## [1] "blue" "red"
##
## [[2]]
## [1] "green" "red"
##
## [[3]]
## [1] "orange" "red"
# simplify = Tを使えば, str_extract_all()が少ないマッチでも最多マッチに合わせた行列を返す
str_extract_all(more, color_match, simplify = TRUE)
## [,1] [,2]
## [1,] "blue" "red"
## [2,] "green" "red"
## [3,] "orange" "red"
x <- c("a", "a b", "a b c")
str_extract_all(x, "[a-z]", simplify = TRUE)
## [,1] [,2] [,3]
## [1,] "a" "" ""
## [2,] "a" "b" ""
## [3,] "a" "b" "c"
色名を単語の途中ではなく、独立した単語として抽出するには、色名を単語の境界を示す\bで挟む。
colors2 <- c("\\bred", "orange", "yellow", "green", "blue", "purple\\b")
color_match2 <- str_c(colors2, collapse = "\\b|\\b")
color_match2
## [1] "\\bred\\b|\\borange\\b|\\byellow\\b|\\bgreen\\b|\\bblue\\b|\\bpurple\\b"
has_color2 <- str_subset(sentences, color_match2)
str_view_all(has_color2, color_match2)
## [1] │ Glue the sheet to the dark <blue> background.
## [2] │ Two <blue> fish swam in the tank.
## [3] │ A wisp of cloud hung in the <blue> air.
## [4] │ Leaves turn brown and <yellow> in the fall.
## [5] │ The spot on the blotter was made by <green> ink.
## [6] │ The sofa cushion is <red> and of light weight.
## [7] │ The sky that morning was clear and bright <blue>.
## [8] │ A <blue> crane is a tall wading bird.
## [9] │ It is hard to erase <blue> or <red> ink.
## [10] │ The lamp shone with a steady <green> flame.
## [11] │ The box is held by a bright <red> snapper.
## [12] │ The houses are built of <red> clay bricks.
## [13] │ The <red> tape bound the smuggled food.
## [14] │ Hedge apples may stain your hands <green>.
## [15] │ The plant grew large and <green> in the window.
## [16] │ The <purple> tie was ten years old.
## [17] │ Bathe and relax in the cool <green> grass.
## [18] │ The lake sparkled in the <red> hot sun.
## [19] │ Mark the spot with a sign painted <red>.
## [20] │ The couch cover and hall drapes were <blue>.
## ... and 9 more
「ing」で終わる単語を抽出する。\\bで単語の区切りを設定。\w*は任意の文字を任意の回数繰り返す。その後にingが続く。
end_ing = "\\b\\w*ing\\b"
# end_ing = "\\b[^ ]+ing\\b" # 別法:"[^ ]+"は「空白以外の文字の1つ以上の連なり」
ing_match <- str_subset(sentences, end_ing)
# str_view_all(ing_match, end_ing)
head(str_extract_all(ing_match, end_ing))
## [[1]]
## [1] "spring"
##
## [[2]]
## [1] "evening"
##
## [[3]]
## [1] "morning"
##
## [[4]]
## [1] "winding"
##
## [[5]]
## [1] "living"
##
## [[6]]
## [1] "king"
文章から名詞を抽出するために、“a”または”the”の後に来る全ての単語を探す。()を使うことにより、2つのマッチがグループになる。
## a または the の後に、空白以外の文字の一つ以上の連なりを名詞の近似とする
noun <- "(a|the) ([^ ]+)"
has_noun <- sentences %>% str_subset(noun) %>% head(10)
has_noun %>% str_extract(noun)
## [1] "the smooth" "the sheet" "the depth" "a chicken" "the parked"
## [6] "the sun" "the huge" "the ball" "the woman" "a helps"
str_extract()は完全マッチ、str_match()は個別の部分マッチを示す、行列を返し第1列が完全マッチ、第2列以降がそれぞれのグループstr_extract()と同様、各文字列全てにマッチしたいなら、str_match_all()が必要。
has_noun %>% str_match(noun)
## [,1] [,2] [,3]
## [1,] "the smooth" "the" "smooth"
## [2,] "the sheet" "the" "sheet"
## [3,] "the depth" "the" "depth"
## [4,] "a chicken" "a" "chicken"
## [5,] "the parked" "the" "parked"
## [6,] "the sun" "the" "sun"
## [7,] "the huge" "the" "huge"
## [8,] "the ball" "the" "ball"
## [9,] "the woman" "the" "woman"
## [10,] "a helps" "a" "helps"
データフレームを検索する場合はtidyr::extractを使用、str_extract()と同様に作用するが、マッチに名前を付ける必要がある。
tibble(sentence = sentences) %>%
tidyr::extract(
sentence, c("article", "noun"), "(a|the) ([^ ]+)",
remove = FALSE
)
## # A tibble: 720 × 3
## sentence article noun
## <chr> <chr> <chr>
## 1 The birch canoe slid on the smooth planks. the smooth
## 2 Glue the sheet to the dark blue background. the sheet
## 3 It's easy to tell the depth of a well. the depth
## 4 These days a chicken leg is a rare dish. a chicken
## 5 Rice is often served in round bowls. <NA> <NA>
## 6 The juice of lemons makes fine punch. <NA> <NA>
## 7 The box was thrown beside the parked truck. the parked
## 8 The hogs were fed chopped corn and garbage. <NA> <NA>
## 9 Four hours of steady work faced us. <NA> <NA>
## 10 A large size in stockings is hard to sell. <NA> <NA>
## # ℹ 710 more rows
x <- c("apple", "pear", "banana")
str_replace(x, "[aeiou]", "-") ## 最初のマッチを置換
## [1] "-pple" "p-ar" "b-nana"
str_replace_all(x, "[aeiou]", "-") ## 全てのマッチを置換
## [1] "-ppl-" "p--r" "b-n-n-"
x <- c("1 house", "2 cars", "3 people")
str_replace_all(x, c("1" = "one", "2" = "two", "3" = "three")) ## 名前付きベクトルを指定すると複数の置換ができる
## [1] "one house" "two cars" "three people"
固定文字列で置換する代わりに、後方参照を使ってマッチした要素を挿入することもできる。
sentences %>%
str_replace("([^ ]+) ([^ ]+) ([^ ]+)", "\\1 \\3 \\2") %>% ## 2番目と3番目の単語の順序を入れ替える
head(5)
## [1] "The canoe birch slid on the smooth planks."
## [2] "Glue sheet the to the dark blue background."
## [3] "It's to easy tell the depth of a well."
## [4] "These a days chicken leg is a rare dish."
## [5] "Rice often is served in round bowls."
str_split()を使うと文字列を分割できる。空白で分割すると、文を単語に分割できる。文ごとに要素の数が異なるのでリストを返す。
sentences %>% head(5) %>% str_split(" ")
## [[1]]
## [1] "The" "birch" "canoe" "slid" "on" "the" "smooth"
## [8] "planks."
##
## [[2]]
## [1] "Glue" "the" "sheet" "to" "the"
## [6] "dark" "blue" "background."
##
## [[3]]
## [1] "It's" "easy" "to" "tell" "the" "depth" "of" "a" "well."
##
## [[4]]
## [1] "These" "days" "a" "chicken" "leg" "is" "a"
## [8] "rare" "dish."
##
## [[5]]
## [1] "Rice" "is" "often" "served" "in" "round" "bowls."
"a|b|c|d" %>%
str_split("\\|") %>% ## 長さ1のベクトルで返したいなら、リストの最初の要素[[1]]だけを抽出
.[[1]]
## [1] "a" "b" "c" "d"
sentences %>% ## simplity =TRUEで行列を返す
head(5) %>%
str_split(" ", simplify = TRUE)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] "The" "birch" "canoe" "slid" "on" "the" "smooth" "planks."
## [2,] "Glue" "the" "sheet" "to" "the" "dark" "blue" "background."
## [3,] "It's" "easy" "to" "tell" "the" "depth" "of" "a"
## [4,] "These" "days" "a" "chicken" "leg" "is" "a" "rare"
## [5,] "Rice" "is" "often" "served" "in" "round" "bowls." ""
## [,9]
## [1,] ""
## [2,] ""
## [3,] "well."
## [4,] "dish."
## [5,] ""
要素の最大個数を指定できる
fields <- c("Name: Hadley", "Country: NZ", "Age: 35")
fields %>% str_split(": ", n = 2, simplify = TRUE)
## [,1] [,2]
## [1,] "Name" "Hadley"
## [2,] "Country" "NZ"
## [3,] "Age" "35"
x <- "This is a sentence. This is another sentence."
str_view_all(x, boundary("word")) ## パターンではなく、文字、行、文、後のboundary()で分割できる
## [1] │ <This> <is> <a> <sentence>. <This> <is> <another> <sentence>.
str_split(x, " ")[[1]]
## [1] "This" "is" "a" "sentence." "" "This"
## [7] "is" "another" "sentence."
str_split(x, boundary("word"))[[1]]
## [1] "This" "is" "a" "sentence" "This" "is" "another"
## [8] "sentence"
文字列のパターンを使えば自動的にregex()にラップされる。
str_view(fruit, "nana")
## [4] │ ba<nana>
# Is shorthand for
str_view(fruit, regex("nana"))
## [4] │ ba<nana>
regex()の他の引数を使ってマッチの詳細を制御できる。
bananas <- c("banana", "Banana", "BANANA")
str_view(bananas, "banana")
## [1] │ <banana>
str_view(bananas, regex("banana", ignore_case = TRUE)) ## ignore_caseで大文字小文字を無視
## [1] │ <banana>
## [2] │ <Banana>
## [3] │ <BANANA>
multilineは^と$が文字列全体の先頭と末尾ではなく、各行の先頭と末尾にマッチする。
x <- "Line 1\nLine 2\nLine 3"
str_extract_all(x, "^Line")[[1]]
## [1] "Line"
str_extract_all(x, regex("^Line", multiline = TRUE))[[1]]
## [1] "Line" "Line" "Line"
commentsは正規表現を読みやすくするためにコメントと空白を使う。
phone <- regex("
\\(? # optional opening parens
(\\d{3}) # area code
[) -]? # optional closing parens, space, or dash
(\\d{3}) # another three numbers
[ -]? # optional space or dash
(\\d{3}) # three more numbers
", comments = TRUE)
str_match("514-791-8141", phone)
## [,1] [,2] [,3] [,4]
## [1,] "514-791-814" "514" "791" "814"
データ解析の対象となるような大きなデータセットは、通常は手作業で入力することは無く、ファイルを読み込んでデータフレームオブジェクトとして保存する。
Environmental paneのImport DatasetボタンからcsvファイルやExcelファイルを読み込むことができる。
read.csv()やread_excel()関数を使用する。
# library(readxl) ## excelやcsv fileを読み込む
# library(openxlsx) ## excelファイルを読み書きする
seiseki <- read.csv("seiseki.csv", header =T, sep = ",", stringsAsFactors = F)
test_data <- read_excel("test_data.xlsx", sheet = "200302") # 読み込むシートも指定できる
head(seiseki)
## kokugo shakai sugaku rika ongaku bijutu taiiku gika eigo
## 1 30 43 51 63 60 66 37 44 20
## 2 39 21 49 56 70 72 56 63 16
## 3 29 30 23 57 69 76 33 54 6
## 4 95 87 77 100 77 82 78 96 87
## 5 70 71 78 67 72 82 46 63 44
## 6 67 53 56 61 61 76 70 66 40
head(test_data)
## # A tibble: 6 × 8
## Date Filename `# of z-slice` defocus `cells with granule`
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 2020-03-02 ade12_006.nd2 22 NA 213
## 2 2020-03-02 ade12_007.nd2 22 NA 258
## 3 2020-03-02 ade12_008.nd2 22 NA 254
## 4 2020-03-02 ade12_020.nd2 22 1 NA
## 5 2020-03-02 ade12_021.nd2 22 NA 307
## 6 2020-03-02 ade12_022.nd2 22 NA 308
## # ℹ 3 more variables: `cell without granule` <dbl>, `Total cell` <dbl>,
## # pct_cell_with_granules <dbl>
Excelファイルやテキストファイル上で該当するデータ範囲をクリップボードにコピーして、データフレームとして読み込みことができる。 WindowsとMacでコマンドが異なる。
# clip_tab <- read.table(pipe("pbpaste"), sep = "", header = T) # Macの場合, sepは区切り文字を指定、header = Tは先頭行を列名にする
# clip_tab <- read.table("clipboard", sep = "", header = T) # Windowsの場合
Rのデータフレームをcsv形式で保存するにはwrite.csv()を、Excel形式で保存するにはwrite.xlsx()をそれぞれ使用。
write.csv(test_data, file = "test_data_csv.csv")
write.xlsx(seiseki, file = "seiseki.xlsx", colNames = T, sheetName = "seiseki")
2つ以上のデータフレームを縦に連結するにはbind_rows()を使用する。同じ名前の列を結合する。
df1 <- read_excel("test_data.xlsx", sheet = "200302")
str(df1)
## tibble [30 × 8] (S3: tbl_df/tbl/data.frame)
## $ Date : chr [1:30] "2020-03-02" "2020-03-02" "2020-03-02" "2020-03-02" ...
## $ Filename : chr [1:30] "ade12_006.nd2" "ade12_007.nd2" "ade12_008.nd2" "ade12_020.nd2" ...
## $ # of z-slice : num [1:30] 22 22 22 22 22 22 22 22 22 22 ...
## $ defocus : num [1:30] NA NA NA 1 NA NA NA NA NA NA ...
## $ cells with granule : num [1:30] 213 258 254 NA 307 308 260 94 144 110 ...
## $ cell without granule : num [1:30] 75 91 109 NA 70 80 55 144 224 181 ...
## $ Total cell : num [1:30] 288 349 363 NA 377 388 315 238 368 291 ...
## $ pct_cell_with_granules: num [1:30] 74 73.9 70 NA 81.4 ...
df2 <- read_excel("test_data.xlsx", sheet = "200304")
str(df2)
## tibble [18 × 8] (S3: tbl_df/tbl/data.frame)
## $ Date : chr [1:18] "2020-03-04" "2020-03-04" "2020-03-04" "2020-03-04" ...
## $ Filename : chr [1:18] "ade12_005.nd2" "ade12_006.nd2" "ade12_007.nd2" "ade12_hyp_008.nd2" ...
## $ # of z-slice : num [1:18] 22 22 22 22 22 22 22 22 22 22 ...
## $ defocus : num [1:18] NA NA NA NA NA NA NA NA NA NA ...
## $ cells with granule : num [1:18] 235 249 192 154 182 146 194 170 120 179 ...
## $ cell without granule : num [1:18] 53 56 57 56 52 52 66 45 50 41 ...
## $ Total cell : num [1:18] 288 305 249 210 234 198 260 215 170 220 ...
## $ pct_cell_with_granules: num [1:18] 81.6 81.6 77.1 73.3 77.8 ...
df_merge <- bind_rows(df1, df2)
str(df_merge)
## tibble [48 × 8] (S3: tbl_df/tbl/data.frame)
## $ Date : chr [1:48] "2020-03-02" "2020-03-02" "2020-03-02" "2020-03-02" ...
## $ Filename : chr [1:48] "ade12_006.nd2" "ade12_007.nd2" "ade12_008.nd2" "ade12_020.nd2" ...
## $ # of z-slice : num [1:48] 22 22 22 22 22 22 22 22 22 22 ...
## $ defocus : num [1:48] NA NA NA 1 NA NA NA NA NA NA ...
## $ cells with granule : num [1:48] 213 258 254 NA 307 308 260 94 144 110 ...
## $ cell without granule : num [1:48] 75 91 109 NA 70 80 55 144 224 181 ...
## $ Total cell : num [1:48] 288 349 363 NA 377 388 315 238 368 291 ...
## $ pct_cell_with_granules: num [1:48] 74 73.9 70 NA 81.4 ...
df_list <- list(df1, df2)
df_merge2 <- bind_rows(df_list) # bind_rowsにわたす引数はリストでも可
多数のデータファイルからデータを読み込んで、1つのデータフレームに集約させる場合は次の3段階で行う:
以下の例ではデータフォルダmaxima_dataの下に5つのサブフォルダがあり、それぞれに複数のcsvファイルが格納されている。
Excelファイルを読み込む場合は.csvを.xlsxに変える。
data_dir <- "maxima_data" # データフォルダの指定
# list.files()関数でフォルダ内のファイル名リストを取得、recursiveでサブフォルダも探索する、拡張子も指定できる
csv_list <- list.files(data_dir, full.names = T, recursive = T, pattern = "*.csv")
str(csv_list) # 合計で84個のcsvファイル
## chr [1:84] "maxima_data/230222_maxima_counts/maxima_counts_Ade4-mNG_sc-lf_003_prom100.csv" ...
apply_res <- csv_list %>%
lapply(., read.csv, header =T, sep = ",", stringsAsFactors = F) # lapplyはリストの各要素に対して演算を行い、結果をリストで返す
length(apply_res) # 84個のデータフレームが格納されている
## [1] 84
csv_merge <- bind_rows(apply_res) # bind_rowsでapply_res中のデータフレームを連結
str(csv_merge)
## 'data.frame': 1260 obs. of 10 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Count : int 29 26 29 39 45 58 69 80 83 101 ...
## $ total_cell_num: int 90 90 90 90 90 90 90 90 90 90 ...
## $ count_per_cell: num 0.322 0.289 0.322 0.433 0.5 ...
## $ prominence : int 100 100 100 100 100 100 100 100 100 100 ...
## $ BC_Mean : num 17.3 17.3 17.3 17.3 17.3 ...
## $ BC_Mode : int 0 0 0 0 0 0 0 0 0 0 ...
## $ BC_Median : int 14 14 14 14 14 14 14 14 14 14 ...
## $ date : int 230222 230222 230222 230222 230222 230222 230222 230222 230222 230222 ...
## $ file : chr "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" ...
# apply_resを作成せずに一気に連結する場合
csv_merge2 <- csv_list %>%
lapply(., read.csv, header =T, sep = ",", stringsAsFactors = F) %>%
bind_rows()
str(csv_merge2)
## 'data.frame': 1260 obs. of 10 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Count : int 29 26 29 39 45 58 69 80 83 101 ...
## $ total_cell_num: int 90 90 90 90 90 90 90 90 90 90 ...
## $ count_per_cell: num 0.322 0.289 0.322 0.433 0.5 ...
## $ prominence : int 100 100 100 100 100 100 100 100 100 100 ...
## $ BC_Mean : num 17.3 17.3 17.3 17.3 17.3 ...
## $ BC_Mode : int 0 0 0 0 0 0 0 0 0 0 ...
## $ BC_Median : int 14 14 14 14 14 14 14 14 14 14 ...
## $ date : int 230222 230222 230222 230222 230222 230222 230222 230222 230222 230222 ...
## $ file : chr "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" "Ade4-mNG_sc-lf_003" ...
読み込んだ多数のデータファイル(.csv or .xlsx)をデータフレームに変換して、1つのExcelファイルのワークシートにそれぞれ出力する場合。
data_dir2 <- "maxima_data/230228_maxima_counts"
csv_list2 <- list.files(data_dir2, full.names = T, recursive = T, pattern = "*.csv")
sheetname_list <- basename(csv_list2) %>% # basenameは上位のフォルダ名を含まないファイル名を返す
str_replace_all(c("maxima_counts_" = "", ".csv" = "")) %>% # str_replace_all()で複数マッチを一度に置換、接頭語と拡張子を除く
str_sub (., start = 1, end = 31) # シート名は31文字以内の制限があるので、さらに先頭から31文字を切り出す
str(sheetname_list)
## chr [1:24] "Ade4-mNG_sc-lf_003_prom100" "Ade4-mNG_sc-lf_003_prom110" ...
names(csv_list2) <- sheetname_list # csv_list2の要素の名前を、短縮したファイル名のリストsheetname_listで置換
file_list <- csv_list2 %>%
lapply(., read.csv, header =T, sep = ",", stringsAsFactors = F) # bind_rowsで結合せずにリストにデータフレームを格納
length(file_list)
## [1] 24
file_list %>% write.xlsx(., file = "assembled.xlsx", colNames = T) ## xlsxファイルに書き出し、各シートが各ファイルに対応
整理データ tidy dataは以下の条件を満たすデータ:
パイプ演算子 %>%
は左辺の計算結果のオブジェクトを右辺の第1引数に渡す。第1引数以外で呼び出す時は
.(ドット演算子)で行う。
# library(magrittr)
3 %>% c(., 3 + .)
## [1] 3 6
3 %>% c(3 + .) # 右辺の第一引数は省略可能
## [1] 3 6
3 %>% { c(3 + .) } # 中括弧はパイプが関数内の最初の引数を使用しないようにする
## [1] 6
mtcars %>% {cor(.$mpg, .$cyl)}
## [1] -0.852162
1:10 %>% mean(.)
## [1] 5.5
1:10 %>% mean # 右辺の第一引数は省略可能
## [1] 5.5
1:10 %>% plot(x = seq(10, 100, 10), y = .)
Dollar演算子 %$%
はデータフレームの列にアクセスできる
iris %>%
filter(Species != "setosa") %$%
lm(Sepal.Length ~ Sepal.Width) %>% ## lm(.$Sepal.Length ~ .$Sepal.Width)と同等
summary()
##
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0032 -0.3877 -0.0774 0.3200 1.7381
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0934 0.4844 6.387 5.70e-09 ***
## Sepal.Width 1.1033 0.1675 6.585 2.27e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5547 on 98 degrees of freedom
## Multiple R-squared: 0.3068, Adjusted R-squared: 0.2997
## F-statistic: 43.36 on 1 and 98 DF, p-value: 2.27e-09
一般的なデータセットは整理データの形式になっていないことの方が多い。
データ解析の前にデータを整理する必要がある。まず何が変数で、何が観測かを見極めることが重要。
日付とその日の天気のような、2つの情報の組を考える。個々の組を識別する情報をKey「キー」と呼び、もう1つをValue「値」と呼ぶ。 Key-Value型データベースは、このようなkeyとvalueの組を保存するためのデータベース、辞書データやハッシュテーブルとしても知られる。 複数のkeyでvalueを識別する場合もある。
列方向すなわち横方向に伸びるようにまとめられたデータをwide format(ワイド形式)、 行方向すなわち縦方向に表が伸びていくようにまとめられたデータをlong format(ロング形式)という。ワイドかロングかは相対的な関係である。 ワイド形式は多数のkeyが列名になっており、valueを示す列は存在しないことが多い。ロング形式は少数のkeyと1つのvalueの列から構成されることが多い。
library(tidyverse)
table4a # tidyverseに収録されているサンプルデータ、ワイド形式、変数の「値」が列名になっている
## # A tibble: 3 × 3
## country `1999` `2000`
## <chr> <dbl> <dbl>
## 1 Afghanistan 745 2666
## 2 Brazil 37737 80488
## 3 China 212258 213766
table4a %>% gather("1999":"2000", key = "year", value = "cases") # gatherでロング形式に変換、変換する列、列名にする変数名key、列名にする値名valueをそれぞれ指定
## # A tibble: 6 × 3
## country year cases
## <chr> <chr> <dbl>
## 1 Afghanistan 1999 745
## 2 Brazil 1999 37737
## 3 China 1999 212258
## 4 Afghanistan 2000 2666
## 5 Brazil 2000 80488
## 6 China 2000 213766
table4a %>% pivot_longer(cols = c("1999", "2000"), names_to = "year", values_to = "cases") # 同機能のより新しい関数
## # A tibble: 6 × 3
## country year cases
## <chr> <chr> <dbl>
## 1 Afghanistan 1999 745
## 2 Afghanistan 2000 2666
## 3 Brazil 1999 37737
## 4 Brazil 2000 80488
## 5 China 1999 212258
## 6 China 2000 213766
table2 # 列typeに変数casesとpopulationが入っており、変数の値は列countに入っている
## # A tibble: 12 × 4
## country year type count
## <chr> <dbl> <chr> <dbl>
## 1 Afghanistan 1999 cases 745
## 2 Afghanistan 1999 population 19987071
## 3 Afghanistan 2000 cases 2666
## 4 Afghanistan 2000 population 20595360
## 5 Brazil 1999 cases 37737
## 6 Brazil 1999 population 172006362
## 7 Brazil 2000 cases 80488
## 8 Brazil 2000 population 174504898
## 9 China 1999 cases 212258
## 10 China 1999 population 1272915272
## 11 China 2000 cases 213766
## 12 China 2000 population 1280428583
table2 %>% spread(key = type, value = count) # spreadで広げる、慣例として既存の列は""で囲まず、新規に作成する際は""で囲む
## # A tibble: 6 × 4
## country year cases population
## <chr> <dbl> <dbl> <dbl>
## 1 Afghanistan 1999 745 19987071
## 2 Afghanistan 2000 2666 20595360
## 3 Brazil 1999 37737 172006362
## 4 Brazil 2000 80488 174504898
## 5 China 1999 212258 1272915272
## 6 China 2000 213766 1280428583
table2 %>% pivot_wider(names_from = type, values_from = count) # 同機能のより新しい関数
## # A tibble: 6 × 4
## country year cases population
## <chr> <dbl> <dbl> <dbl>
## 1 Afghanistan 1999 745 19987071
## 2 Afghanistan 2000 2666 20595360
## 3 Brazil 1999 37737 172006362
## 4 Brazil 2000 80488 174504898
## 5 China 1999 212258 1272915272
## 6 China 2000 213766 1280428583
1つの列に2つの変数が含まれている場合、separate()で分ける
table3
## # A tibble: 6 × 3
## country year rate
## <chr> <dbl> <chr>
## 1 Afghanistan 1999 745/19987071
## 2 Afghanistan 2000 2666/20595360
## 3 Brazil 1999 37737/172006362
## 4 Brazil 2000 80488/174504898
## 5 China 1999 212258/1272915272
## 6 China 2000 213766/1280428583
table3 %>% separate(rate, into = c("cases", "population")) ## separate()はdefaultで非英数文字で値を分割する
## # A tibble: 6 × 4
## country year cases population
## <chr> <dbl> <chr> <chr>
## 1 Afghanistan 1999 745 19987071
## 2 Afghanistan 2000 2666 20595360
## 3 Brazil 1999 37737 172006362
## 4 Brazil 2000 80488 174504898
## 5 China 1999 212258 1272915272
## 6 China 2000 213766 1280428583
table3 %>% separate(rate, into = c("cases", "population"), sep = "/") ## sep引数で指定することも可能
## # A tibble: 6 × 4
## country year cases population
## <chr> <dbl> <chr> <chr>
## 1 Afghanistan 1999 745 19987071
## 2 Afghanistan 2000 2666 20595360
## 3 Brazil 1999 37737 172006362
## 4 Brazil 2000 80488 174504898
## 5 China 1999 212258 1272915272
## 6 China 2000 213766 1280428583
table3 %>% separate(rate, into = c("cases", "population"), convert = T) ## convertにより数字がより良いint型に変換される
## # A tibble: 6 × 4
## country year cases population
## <chr> <dbl> <int> <int>
## 1 Afghanistan 1999 745 19987071
## 2 Afghanistan 2000 2666 20595360
## 3 Brazil 1999 37737 172006362
## 4 Brazil 2000 80488 174504898
## 5 China 1999 212258 1272915272
## 6 China 2000 213766 1280428583
table3 %>% separate(year, into = c("century", "year"), sep = 2) ## sepに整数値を渡すと、文字列の左端を+1, 右端を-1としてその整数の位置で分割する
## # A tibble: 6 × 4
## country century year rate
## <chr> <chr> <chr> <chr>
## 1 Afghanistan 19 99 745/19987071
## 2 Afghanistan 20 00 2666/20595360
## 3 Brazil 19 99 37737/172006362
## 4 Brazil 20 00 80488/174504898
## 5 China 19 99 212258/1272915272
## 6 China 20 00 213766/1280428583
unite()はseparate()の逆の変換。
table5
## # A tibble: 6 × 4
## country century year rate
## <chr> <chr> <chr> <chr>
## 1 Afghanistan 19 99 745/19987071
## 2 Afghanistan 20 00 2666/20595360
## 3 Brazil 19 99 37737/172006362
## 4 Brazil 20 00 80488/174504898
## 5 China 19 99 212258/1272915272
## 6 China 20 00 213766/1280428583
table5 %>% unite("new", century, year, sep = "") ## sepを指定しないデフォルトでは「_」で値が繋がれる
## # A tibble: 6 × 3
## country new rate
## <chr> <chr> <chr>
## 1 Afghanistan 1999 745/19987071
## 2 Afghanistan 2000 2666/20595360
## 3 Brazil 1999 37737/172006362
## 4 Brazil 2000 80488/174504898
## 5 China 1999 212258/1272915272
## 6 China 2000 213766/1280428583
dplyrの基本関数:
filter():値から観測値を選び出すarrange():行を並び替えるselect():名前で変数(列)を選ぶmutate():既存の変数の関数で新たな変数を作るsummarise():多数の値から単一の要約量を作るlibrary(tidyverse)
library(nycflights13)
data("flights")
filter(flights, month==1, day==1) ## 1月1日のフライト、複数の条件は&になるので両方満たした行が取り出される
## # A tibble: 842 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 -1 1004 1022
## 5 2013 1 1 554 600 -6 812 837
## 6 2013 1 1 554 558 -4 740 728
## 7 2013 1 1 555 600 -5 913 854
## 8 2013 1 1 557 600 -3 709 723
## 9 2013 1 1 557 600 -3 838 846
## 10 2013 1 1 558 600 -2 753 745
## # ℹ 832 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
filter(flights, month == 11 | month ==12) ## 11月または12月のフライト
## # A tibble: 55,403 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 11 1 5 2359 6 352 345
## 2 2013 11 1 35 2250 105 123 2356
## 3 2013 11 1 455 500 -5 641 651
## 4 2013 11 1 539 545 -6 856 827
## 5 2013 11 1 542 545 -3 831 855
## 6 2013 11 1 549 600 -11 912 923
## 7 2013 11 1 550 600 -10 705 659
## 8 2013 11 1 554 600 -6 659 701
## 9 2013 11 1 554 600 -6 826 827
## 10 2013 11 1 554 600 -6 749 751
## # ℹ 55,393 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
# nov_dec <- filter(flights, month %in% c(11, 12)) ## このようにも書ける
filter(flights, arr_delay <= 120, dep_delay <= 120) ## 出発または到着で120分を超える遅延がなかったフライト
## # A tibble: 316,050 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 -1 1004 1022
## 5 2013 1 1 554 600 -6 812 837
## 6 2013 1 1 554 558 -4 740 728
## 7 2013 1 1 555 600 -5 913 854
## 8 2013 1 1 557 600 -3 709 723
## 9 2013 1 1 557 600 -3 838 846
## 10 2013 1 1 558 600 -2 753 745
## # ℹ 316,040 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
# filter(flights, !(arr_delay > 120 | dep_delay > 120)) # このようにも書ける
filter(flights, between(dep_time, 0, 6)) ## betweenは範囲を指定、指定した数値も含む
## # A tibble: 155 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 9 2 2359 3 432 444
## 2 2013 1 10 3 2359 4 426 437
## 3 2013 1 13 1 2249 72 108 2357
## 4 2013 1 13 2 2359 3 502 444
## 5 2013 1 13 3 2030 213 340 2350
## 6 2013 1 16 2 2125 157 119 2250
## 7 2013 1 16 3 1946 257 212 2154
## 8 2013 1 22 5 2359 6 452 444
## 9 2013 1 30 3 2159 124 100 2306
## 10 2013 1 31 1 2100 181 124 2225
## # ℹ 145 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
# 順序付ける列名の集合を引数に取る。
# 複数の列名を与えると、前列の順序で同じ順序になったものをさらに順序付けるために後列を使用。
arrange(flights, year, month, day)
## # A tibble: 336,776 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 -1 1004 1022
## 5 2013 1 1 554 600 -6 812 837
## 6 2013 1 1 554 558 -4 740 728
## 7 2013 1 1 555 600 -5 913 854
## 8 2013 1 1 557 600 -3 709 723
## 9 2013 1 1 557 600 -3 838 846
## 10 2013 1 1 558 600 -2 753 745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
arrange(flights, desc(arr_delay)) ## 降順にするにはdescを使用、NAは常に一番最後になる
## # A tibble: 336,776 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 9 641 900 1301 1242 1530
## 2 2013 6 15 1432 1935 1137 1607 2120
## 3 2013 1 10 1121 1635 1126 1239 1810
## 4 2013 9 20 1139 1845 1014 1457 2210
## 5 2013 7 22 845 1600 1005 1044 1815
## 6 2013 4 10 1100 1900 960 1342 2211
## 7 2013 3 17 2321 810 911 135 1020
## 8 2013 7 22 2257 759 898 121 1026
## 9 2013 12 5 756 1700 896 1058 2020
## 10 2013 5 3 1133 2055 878 1250 2215
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
select(flights, year, month, day) ## 列を名前で選ぶ
## # A tibble: 336,776 × 3
## year month day
## <int> <int> <int>
## 1 2013 1 1
## 2 2013 1 1
## 3 2013 1 1
## 4 2013 1 1
## 5 2013 1 1
## 6 2013 1 1
## 7 2013 1 1
## 8 2013 1 1
## 9 2013 1 1
## 10 2013 1 1
## # ℹ 336,766 more rows
select(flights, year:day) ## yearとdayの間にある全ての列を選ぶ
## # A tibble: 336,776 × 3
## year month day
## <int> <int> <int>
## 1 2013 1 1
## 2 2013 1 1
## 3 2013 1 1
## 4 2013 1 1
## 5 2013 1 1
## 6 2013 1 1
## 7 2013 1 1
## 8 2013 1 1
## 9 2013 1 1
## 10 2013 1 1
## # ℹ 336,766 more rows
select(flights, -(year:day)) ## yearとdayの間の列以外の全ての列
## # A tibble: 336,776 × 16
## dep_time sched_dep_time dep_delay arr_time sched_arr_time arr_delay carrier
## <int> <int> <dbl> <int> <int> <dbl> <chr>
## 1 517 515 2 830 819 11 UA
## 2 533 529 4 850 830 20 UA
## 3 542 540 2 923 850 33 AA
## 4 544 545 -1 1004 1022 -18 B6
## 5 554 600 -6 812 837 -25 DL
## 6 554 558 -4 740 728 12 UA
## 7 555 600 -5 913 854 19 B6
## 8 557 600 -3 709 723 -14 EV
## 9 557 600 -3 838 846 -8 B6
## 10 558 600 -2 753 745 8 AA
## # ℹ 336,766 more rows
## # ℹ 9 more variables: flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## # air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>
select()ではヘルパー関数
start_with()、ends_with()、contains()、num_range()等が使用できる。
select(flights, starts_with("dep")) # "dep"で始まる列
## # A tibble: 336,776 × 2
## dep_time dep_delay
## <int> <dbl>
## 1 517 2
## 2 533 4
## 3 542 2
## 4 544 -1
## 5 554 -6
## 6 554 -4
## 7 555 -5
## 8 557 -3
## 9 557 -3
## 10 558 -2
## # ℹ 336,766 more rows
select(flights, ends_with("time")) # "time"で終わる列
## # A tibble: 336,776 × 5
## dep_time sched_dep_time arr_time sched_arr_time air_time
## <int> <int> <int> <int> <dbl>
## 1 517 515 830 819 227
## 2 533 529 850 830 227
## 3 542 540 923 850 160
## 4 544 545 1004 1022 183
## 5 554 600 812 837 116
## 6 554 558 740 728 150
## 7 555 600 913 854 158
## 8 557 600 709 723 53
## 9 557 600 838 846 140
## 10 558 600 753 745 138
## # ℹ 336,766 more rows
select(flights, contains(c("air", "sched", "arr"))) # "air", "sched", "arr"のいずれかを含む
## # A tibble: 336,776 × 6
## air_time sched_dep_time sched_arr_time arr_time arr_delay carrier
## <dbl> <int> <int> <int> <dbl> <chr>
## 1 227 515 819 830 11 UA
## 2 227 529 830 850 20 UA
## 3 160 540 850 923 33 AA
## 4 183 545 1022 1004 -18 B6
## 5 116 600 837 812 -25 DL
## 6 150 558 728 740 12 UA
## 7 158 600 854 913 19 B6
## 8 53 600 723 709 -14 EV
## 9 140 600 846 838 -8 B6
## 10 138 600 745 753 8 AA
## # ℹ 336,766 more rows
select(flights, matches("time$")) # 正規表現にマッチする列名、ここでは"time"で終わる列
## # A tibble: 336,776 × 5
## dep_time sched_dep_time arr_time sched_arr_time air_time
## <int> <int> <int> <int> <dbl>
## 1 517 515 830 819 227
## 2 533 529 850 830 227
## 3 542 540 923 850 160
## 4 544 545 1004 1022 183
## 5 554 600 812 837 116
## 6 554 558 740 728 150
## 7 555 600 913 854 158
## 8 557 600 709 723 53
## 9 557 600 838 846 140
## 10 558 600 753 745 138
## # ℹ 336,766 more rows
# num_range("x", 1:3) ## x1, x2, x3にマッチする
select()で列名を変更することもできるが、明示しない変数を全て削除するのであまり便利でない。
列名の変更にはrename()を使用する。
select(flights, tail_num = tailnum) # new_name = old_name
## # A tibble: 336,776 × 1
## tail_num
## <chr>
## 1 N14228
## 2 N24211
## 3 N619AA
## 4 N804JB
## 5 N668DN
## 6 N39463
## 7 N516JB
## 8 N829AS
## 9 N593JB
## 10 N3ALAA
## # ℹ 336,766 more rows
rename(flights, tail_num = tailnum)
## # A tibble: 336,776 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 -1 1004 1022
## 5 2013 1 1 554 600 -6 812 837
## 6 2013 1 1 554 558 -4 740 728
## 7 2013 1 1 555 600 -5 913 854
## 8 2013 1 1 557 600 -3 709 723
## 9 2013 1 1 557 600 -3 838 846
## 10 2013 1 1 558 600 -2 753 745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tail_num <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
select()とヘルパー関数everything()を同時に使うことで、指定した列を先頭に移動できる。
select(flights, time_hour, air_time, everything())
## # A tibble: 336,776 × 19
## time_hour air_time year month day dep_time sched_dep_time
## <dttm> <dbl> <int> <int> <int> <int> <int>
## 1 2013-01-01 05:00:00 227 2013 1 1 517 515
## 2 2013-01-01 05:00:00 227 2013 1 1 533 529
## 3 2013-01-01 05:00:00 160 2013 1 1 542 540
## 4 2013-01-01 05:00:00 183 2013 1 1 544 545
## 5 2013-01-01 06:00:00 116 2013 1 1 554 600
## 6 2013-01-01 05:00:00 150 2013 1 1 554 558
## 7 2013-01-01 06:00:00 158 2013 1 1 555 600
## 8 2013-01-01 06:00:00 53 2013 1 1 557 600
## 9 2013-01-01 06:00:00 140 2013 1 1 557 600
## 10 2013-01-01 06:00:00 138 2013 1 1 558 600
## # ℹ 336,766 more rows
## # ℹ 12 more variables: dep_delay <dbl>, arr_time <int>, sched_arr_time <int>,
## # arr_delay <dbl>, carrier <chr>, flight <int>, tailnum <chr>, origin <chr>,
## # dest <chr>, distance <dbl>, hour <dbl>, minute <dbl>
選択したい列名を文字列ベクトルで指定したい場合はone_of()を使用する。
vars <- c("year", "month", "day", "dep_delay", "arr_delay")
select(flights, one_of(vars)) ## one_of()は文字列をselectにわたす場合に必要
## # A tibble: 336,776 × 5
## year month day dep_delay arr_delay
## <int> <int> <int> <dbl> <dbl>
## 1 2013 1 1 2 11
## 2 2013 1 1 4 20
## 3 2013 1 1 2 33
## 4 2013 1 1 -1 -18
## 5 2013 1 1 -6 -25
## 6 2013 1 1 -4 12
## 7 2013 1 1 -5 19
## 8 2013 1 1 -3 -14
## 9 2013 1 1 -3 -8
## 10 2013 1 1 -2 8
## # ℹ 336,766 more rows
# select(flights, vars) # 直接使用すると警告が出る
flights_sml <- select(flights, ## 変数が少ないデータセットを作る
year:day, ends_with("delay"),
distance,
air_time)
mutate(flights_sml,
gain = arr_delay - dep_delay, # 新しく定義した列を最後尾に追加する
speed = distance/air_time *60,
hours = air_time*60,
gain_per_hour = gain/hours) ## 作成したばかりの列を参照できる
## # A tibble: 336,776 × 11
## year month day dep_delay arr_delay distance air_time gain speed hours
## <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2013 1 1 2 11 1400 227 9 370. 13620
## 2 2013 1 1 4 20 1416 227 16 374. 13620
## 3 2013 1 1 2 33 1089 160 31 408. 9600
## 4 2013 1 1 -1 -18 1576 183 -17 517. 10980
## 5 2013 1 1 -6 -25 762 116 -19 394. 6960
## 6 2013 1 1 -4 12 719 150 16 288. 9000
## 7 2013 1 1 -5 19 1065 158 24 404. 9480
## 8 2013 1 1 -3 -14 229 53 -11 259. 3180
## 9 2013 1 1 -3 -8 944 140 -5 405. 8400
## 10 2013 1 1 -2 8 733 138 10 319. 8280
## # ℹ 336,766 more rows
## # ℹ 1 more variable: gain_per_hour <dbl>
普通の書き方では新しく作る列名には変数を使用できず、変数名がそのまま新規の列名として解釈される。
flights_sml2 <- select(flights, ## 変数が少ないデータセットを作る
month,
distance,
air_time)
for (i in 1:3){
new_col_name = paste("new_col", i, sep = "_")
flights_sml2 <-
flights_sml2 %>%
mutate(new_col_name = n())
}
flights_sml2
## # A tibble: 336,776 × 4
## month distance air_time new_col_name
## <int> <dbl> <dbl> <int>
## 1 1 1400 227 336776
## 2 1 1416 227 336776
## 3 1 1089 160 336776
## 4 1 1576 183 336776
## 5 1 762 116 336776
## 6 1 719 150 336776
## 7 1 1065 158 336776
## 8 1 229 53 336776
## 9 1 944 140 336776
## 10 1 733 138 336776
## # ℹ 336,766 more rows
列名に変数を使用するには先頭に!!(読み方:bang
bang)を付けると変数として解釈される。またこの場合は代入に=ではなく:=を使用する必要がある。
for (i in 1:3){
new_col_name = paste("new_col", i, sep = "_")
flights_sml2 <-
flights_sml2 %>%
mutate(!!new_col_name := n())
}
flights_sml2
## # A tibble: 336,776 × 7
## month distance air_time new_col_name new_col_1 new_col_2 new_col_3
## <int> <dbl> <dbl> <int> <int> <int> <int>
## 1 1 1400 227 336776 336776 336776 336776
## 2 1 1416 227 336776 336776 336776 336776
## 3 1 1089 160 336776 336776 336776 336776
## 4 1 1576 183 336776 336776 336776 336776
## 5 1 762 116 336776 336776 336776 336776
## 6 1 719 150 336776 336776 336776 336776
## 7 1 1065 158 336776 336776 336776 336776
## 8 1 229 53 336776 336776 336776 336776
## 9 1 944 140 336776 336776 336776 336776
## 10 1 733 138 336776 336776 336776 336776
## # ℹ 336,766 more rows
transmute()は作成した列を抽出する
transmute(flights,
gain = arr_delay - dep_delay,
hours = 60*air_time,
gain_per_hours = gain/hours)
## # A tibble: 336,776 × 3
## gain hours gain_per_hours
## <dbl> <dbl> <dbl>
## 1 9 13620 0.000661
## 2 16 13620 0.00117
## 3 31 9600 0.00323
## 4 -17 10980 -0.00155
## 5 -19 6960 -0.00273
## 6 16 9000 0.00178
## 7 24 9480 0.00253
## 8 -11 3180 -0.00346
## 9 -5 8400 -0.000595
## 10 10 8280 0.00121
## # ℹ 336,766 more rows
summarise()はデータフレームを要約して結果を1行で出力する。後述のgroup_by()と組み合わせて使用することが多い。
便利な要約関数
median()sd(), IOR(), mad()min(), max(), quantile()first(), nth(x, 2), last()summarise(flights, mean_delay= mean(dep_delay, na.rm =T), max_delay = max(dep_delay, na.rm = T)) ## 欠損値は除く
## # A tibble: 1 × 2
## mean_delay max_delay
## <dbl> <dbl>
## 1 12.6 1301
group_by()はデータを指定されたカテゴリ変数の値ごとにグループ分けし、分析単位をデータセット全体から個別グループに変更する。
flights %>% group_by(year, month) %>%
summarise(mean_dep_delay=mean(dep_delay, na.rm =T)) # 各年の各月ごとのdep_delayの平均値を計算
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 3
## # Groups: year [1]
## year month mean_dep_delay
## <int> <int> <dbl>
## 1 2013 1 10.0
## 2 2013 2 10.8
## 3 2013 3 13.2
## 4 2013 4 13.9
## 5 2013 5 13.0
## 6 2013 6 20.8
## 7 2013 7 21.7
## 8 2013 8 12.6
## 9 2013 9 6.72
## 10 2013 10 6.24
## 11 2013 11 5.44
## 12 2013 12 16.6
flights %>% group_by(carrier) %>%
filter(dep_delay >= 0) %>% # グループごとにfilter()を適用
summarise(mean_dep_delay=mean(dep_delay, na.rm =T))
## # A tibble: 16 × 2
## carrier mean_dep_delay
## <chr> <dbl>
## 1 9E 44.9
## 2 AA 32.1
## 3 AS 27.9
## 4 B6 35.2
## 5 DL 31.5
## 6 EV 47.0
## 7 F9 40.0
## 8 FL 37.8
## 9 HA 37.3
## 10 MQ 38.6
## 11 OO 58
## 12 UA 26.6
## 13 US 29.2
## 14 VX 29.1
## 15 WN 30.3
## 16 YV 49.2
flights %>%
group_by(dest) %>%
summarise(
count = n(),
dist = mean(distance, na.rm=T),
delay = mean(arr_delay, na.rm=T)
) %>%
filter(count >20, dest != "HML")
## # A tibble: 97 × 4
## dest count dist delay
## <chr> <int> <dbl> <dbl>
## 1 ABQ 254 1826 4.38
## 2 ACK 265 199 4.85
## 3 ALB 439 143 14.4
## 4 ATL 17215 757. 11.3
## 5 AUS 2439 1514. 6.02
## 6 AVL 275 584. 8.00
## 7 BDL 443 116 7.05
## 8 BGR 375 378 8.03
## 9 BHM 297 866. 16.9
## 10 BNA 6333 758. 11.8
## # ℹ 87 more rows
n()は各グループの行数(サイズ)を返す。非欠損値の個数を数えるにはsum(!is.na(x))を使用し、一意な値の個数を数えるにはn_distinct()を使う。
flights %>%
group_by(carrier, origin) %>%
summarise(count = n(), unique_count = n_distinct(tailnum))
## `summarise()` has grouped output by 'carrier'. You can override using the
## `.groups` argument.
## # A tibble: 35 × 4
## # Groups: carrier [16]
## carrier origin count unique_count
## <chr> <chr> <int> <int>
## 1 9E EWR 1268 199
## 2 9E JFK 14651 204
## 3 9E LGA 2541 192
## 4 AA EWR 3487 481
## 5 AA JFK 13783 409
## 6 AA LGA 15459 431
## 7 AS EWR 714 84
## 8 B6 EWR 6557 190
## 9 B6 JFK 42076 193
## 10 B6 LGA 6002 128
## # ℹ 25 more rows
count(a, b)はdf %>% group_by(a, b) %>% summarise(count = n())に等しい。
オプションで重み変数wtを追加することもできる。
flights %>% count(year, dest)
## # A tibble: 105 × 3
## year dest n
## <int> <chr> <int>
## 1 2013 ABQ 254
## 2 2013 ACK 265
## 3 2013 ALB 439
## 4 2013 ANC 8
## 5 2013 ATL 17215
## 6 2013 AUS 2439
## 7 2013 AVL 275
## 8 2013 BDL 443
## 9 2013 BGR 375
## 10 2013 BHM 297
## # ℹ 95 more rows
flights %>% count(tailnum, wt = distance) # 個別の飛行機の総飛行距離を「カウント」できる
## # A tibble: 4,044 × 2
## tailnum n
## <chr> <dbl>
## 1 D942DN 3418
## 2 N0EGMQ 250866
## 3 N10156 115966
## 4 N102UW 25722
## 5 N103US 24619
## 6 N104UW 25157
## 7 N10575 150194
## 8 N105UW 23618
## 9 N107US 21677
## 10 N108UW 32070
## # ℹ 4,034 more rows
最頻値を直接求める関数は無いので、個数や件数をカウントしてからその最大値を抽出する。
例:月ごとの目的地(dest)の最頻値を求める。
monthとdestでグループ化して個数countを求める。さらにcountの最大値max_countを求め、それとcountが一致する行を抽出。
flights %>%
group_by(month, dest) %>%
summarise(count = n()) %>%
mutate(max_count = max(count)) %>%
filter(count == max_count)
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 4
## # Groups: month [12]
## month dest count max_count
## <int> <chr> <int> <int>
## 1 1 ATL 1396 1396
## 2 2 ATL 1267 1267
## 3 3 ATL 1448 1448
## 4 4 ATL 1490 1490
## 5 5 ORD 1582 1582
## 6 6 ORD 1547 1547
## 7 7 ORD 1573 1573
## 8 8 ORD 1604 1604
## 9 9 ORD 1582 1582
## 10 10 ORD 1604 1604
## 11 11 ATL 1384 1384
## 12 12 ATL 1463 1463
countの最大値はcount %>% max()と書くこともできる。
flights %>%
group_by(month, dest) %>%
summarise(count = n()) %>%
filter(count == count %>% max())
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 3
## # Groups: month [12]
## month dest count
## <int> <chr> <int>
## 1 1 ATL 1396
## 2 2 ATL 1267
## 3 3 ATL 1448
## 4 4 ATL 1490
## 5 5 ORD 1582
## 6 6 ORD 1547
## 7 7 ORD 1573
## 8 8 ORD 1604
## 9 9 ORD 1582
## 10 10 ORD 1604
## 11 11 ATL 1384
## 12 12 ATL 1463
クロス集計表を使用する方法。
table()でクロス集計が可能。
table_dest <- table(flights$month, flights$dest)
table_dest[, 1:10]
##
## ABQ ACK ALB ANC ATL AUS AVL BDL BGR BHM
## 1 0 0 64 0 1396 169 2 37 0 25
## 2 0 0 58 0 1267 167 0 46 0 23
## 3 0 0 57 0 1448 273 0 62 2 26
## 4 9 0 13 0 1490 209 3 61 29 26
## 5 31 21 59 0 1499 204 35 57 29 27
## 6 30 43 34 0 1438 203 40 48 28 27
## 7 31 66 15 4 1511 213 40 2 30 29
## 8 31 67 20 4 1507 211 35 50 30 30
## 9 30 45 20 0 1364 205 30 39 56 25
## 10 31 23 1 0 1448 214 31 9 60 27
## 11 30 0 46 0 1384 186 30 3 55 22
## 12 31 0 52 0 1463 185 29 29 56 10
forループでクロス集計表のi行目の最大の要素の目的地(dest)を取り出す。which.max()は最大値のインデックスを返す。
months <- names(table_dest[,1]) # 1列目を取り出す
mode_dest <- c()
for (i in 1:length(months)){
mode_dest[i] <- names(which.max(table_dest[i,])) # i行目のデータで最大の要素の名前を取り出す
}
data.frame(month = months, mode_dest = mode_dest)
## month mode_dest
## 1 1 ATL
## 2 2 ATL
## 3 3 ATL
## 4 4 ATL
## 5 5 ORD
## 6 6 ORD
## 7 7 ORD
## 8 8 ORD
## 9 9 ORD
## 10 10 ORD
## 11 11 ATL
## 12 12 ATL
filter(), arrange(), select(), mutate(), summarise(), group_by()には複数の列を同時に処理する、scoped
functionが用意されている。
_all: 全ての列に同じ操作を行う。_if: 条件を満たす列に同じ操作を行う。_at: 指定した複数の列に同じ操作を行う。flights %>% select(ends_with("time")) %>% # "time"で終わる列を選択
mutate_all(abs) # 関数abs()で全ての列の絶対値を求める
## # A tibble: 336,776 × 5
## dep_time sched_dep_time arr_time sched_arr_time air_time
## <int> <int> <int> <int> <dbl>
## 1 517 515 830 819 227
## 2 533 529 850 830 227
## 3 542 540 923 850 160
## 4 544 545 1004 1022 183
## 5 554 600 812 837 116
## 6 554 558 740 728 150
## 7 555 600 913 854 158
## 8 557 600 709 723 53
## 9 557 600 838 846 140
## 10 558 600 753 745 138
## # ℹ 336,766 more rows
flights %>% mutate_if(is.numeric, abs) # 数値型のデータ列に対してのみ、abs()を適用
## # A tibble: 336,776 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 1 1004 1022
## 5 2013 1 1 554 600 6 812 837
## 6 2013 1 1 554 558 4 740 728
## 7 2013 1 1 555 600 5 913 854
## 8 2013 1 1 557 600 3 709 723
## 9 2013 1 1 557 600 3 838 846
## 10 2013 1 1 558 600 2 753 745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
flights %>% mutate_at(c("dep_delay", "arr_delay"), abs) # 列を直接指定できる
## # A tibble: 336,776 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 1 1004 1022
## 5 2013 1 1 554 600 6 812 837
## 6 2013 1 1 554 558 4 740 728
## 7 2013 1 1 555 600 5 913 854
## 8 2013 1 1 557 600 3 709 723
## 9 2013 1 1 557 600 3 838 846
## 10 2013 1 1 558 600 2 753 745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
最近ではscoped
functionの代わりにacross()の使用が推奨されている。
flights %>% select(ends_with("time")) %>% # "time"で終わる列を選択
mutate(across(everything(), abs)) # everything()で全ての列を指定
## # A tibble: 336,776 × 5
## dep_time sched_dep_time arr_time sched_arr_time air_time
## <int> <int> <int> <int> <dbl>
## 1 517 515 830 819 227
## 2 533 529 850 830 227
## 3 542 540 923 850 160
## 4 544 545 1004 1022 183
## 5 554 600 812 837 116
## 6 554 558 740 728 150
## 7 555 600 913 854 158
## 8 557 600 709 723 53
## 9 557 600 838 846 140
## 10 558 600 753 745 138
## # ℹ 336,766 more rows
flights %>% mutate(across(is.numeric, abs)) # 数値型のデータ列に対してのみ、abs()を適用
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(is.numeric, abs)`.
## Caused by warning:
## ! Use of bare predicate functions was deprecated in tidyselect 1.1.0.
## ℹ Please use wrap predicates in `where()` instead.
## # Was:
## data %>% select(is.numeric)
##
## # Now:
## data %>% select(where(is.numeric))
## # A tibble: 336,776 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 1 1004 1022
## 5 2013 1 1 554 600 6 812 837
## 6 2013 1 1 554 558 4 740 728
## 7 2013 1 1 555 600 5 913 854
## 8 2013 1 1 557 600 3 709 723
## 9 2013 1 1 557 600 3 838 846
## 10 2013 1 1 558 600 2 753 745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
flights %>% mutate(across(c("dep_delay", "arr_delay"), abs)) # 列を直接指定できる
## # A tibble: 336,776 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 1 1004 1022
## 5 2013 1 1 554 600 6 812 837
## 6 2013 1 1 554 558 4 740 728
## 7 2013 1 1 555 600 5 913 854
## 8 2013 1 1 557 600 3 709 723
## 9 2013 1 1 557 600 3 838 846
## 10 2013 1 1 558 600 2 753 745
## # ℹ 336,766 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
flights %>% summarise(across(ends_with("time"), ~ mean(.x, na.rm = TRUE))) # ラムダ関数も使える
## # A tibble: 1 × 5
## dep_time sched_dep_time arr_time sched_arr_time air_time
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1349. 1344. 1502. 1536. 151.
dplyr::lag()やdplyr::lead()により変数(ベクトル)の前後の値を取得して、差分を簡単に計算することができる。ずらすことにより欠損値になる箇所はdefaultで指定した値で埋められる。
単にlag()と書くとstats::lag()関数になるので注意。
x <- 1:10
lag1 <- dplyr::lag(x, n = 1, default = NA) # nはずらすindexの大きさ、
lead1 <- dplyr::lead(x)
data.frame(x, lag1, lead1)
## x lag1 lead1
## 1 1 NA 2
## 2 2 1 3
## 3 3 2 4
## 4 4 3 5
## 5 5 4 6
## 6 6 5 7
## 7 7 6 8
## 8 8 7 9
## 9 9 8 10
## 10 10 9 NA
月ごとの総飛行距離の前月との差分を計算。
flights %>%
group_by(month) %>%
summarise(sum_distance = sum(distance)) %>%
mutate(lag_month = lag(month), lag_sum_distance = lag(sum_distance),
diff = sum_distance - lag_sum_distance)
## # A tibble: 12 × 5
## month sum_distance lag_month lag_sum_distance diff
## <int> <dbl> <int> <dbl> <dbl>
## 1 1 27188805 NA NA NA
## 2 2 24975509 1 27188805 -2213296
## 3 3 29179636 2 24975509 4204127
## 4 4 29427294 3 29179636 247658
## 5 5 29974128 4 29427294 546834
## 6 6 29856388 5 29974128 -117740
## 7 7 31149199 6 29856388 1292811
## 8 8 31149334 7 31149199 135
## 9 9 28711426 8 31149334 -2437908
## 10 10 30012086 9 28711426 1300660
## 11 11 28639718 10 30012086 -1372368
## 12 12 29954084 11 28639718 1314366
dplyr::distinct()によりデータフレームから重複を削除できる。引数を指定しないと重複の無いデータフレームを返す。
# originとdestの組み合わせ
flights %>%
select(origin, dest) %>%
distinct()
## # A tibble: 224 × 2
## origin dest
## <chr> <chr>
## 1 EWR IAH
## 2 LGA IAH
## 3 JFK MIA
## 4 JFK BQN
## 5 LGA ATL
## 6 EWR ORD
## 7 EWR FLL
## 8 LGA IAD
## 9 JFK MCO
## 10 LGA ORD
## # ℹ 214 more rows
引数に列名を指定することで、重複を確認する列を指定することができる。列名は複数指定することも可能。
flights %>%
select(origin, dest) %>%
distinct(origin)
## # A tibble: 3 × 1
## origin
## <chr>
## 1 EWR
## 2 LGA
## 3 JFK
flights %>%
select(origin, dest) %>%
distinct(dest)
## # A tibble: 105 × 1
## dest
## <chr>
## 1 IAH
## 2 MIA
## 3 BQN
## 4 ATL
## 5 ORD
## 6 FLL
## 7 IAD
## 8 MCO
## 9 PBI
## 10 TPA
## # ℹ 95 more rows
flights %>%
select(origin, dest) %>%
distinct(origin, dest) # distinct()と同じ
## # A tibble: 224 × 2
## origin dest
## <chr> <chr>
## 1 EWR IAH
## 2 LGA IAH
## 3 JFK MIA
## 4 JFK BQN
## 5 LGA ATL
## 6 EWR ORD
## 7 EWR FLL
## 8 LGA IAD
## 9 JFK MCO
## 10 LGA ORD
## # ℹ 214 more rows
引数.keep_all = TRUEとすると指定した列以外の全変数も残す。この時残るのは初出の行。
flights %>%
select(origin, dest) %>%
distinct(origin, .keep_all = TRUE)
## # A tibble: 3 × 2
## origin dest
## <chr> <chr>
## 1 EWR IAH
## 2 LGA IAH
## 3 JFK MIA
関係データベース(Relational
database)とは表形式の複数のデータを関連付けて扱うデータベース。
表の組を結びつける変数はkeyと呼ばれる。keyは1つまたは複数で観測を一意に表し、以下の2種類がある:
library(tidyverse)
library(nycflights13)
airlines
## # A tibble: 16 × 2
## carrier name
## <chr> <chr>
## 1 9E Endeavor Air Inc.
## 2 AA American Airlines Inc.
## 3 AS Alaska Airlines Inc.
## 4 B6 JetBlue Airways
## 5 DL Delta Air Lines Inc.
## 6 EV ExpressJet Airlines Inc.
## 7 F9 Frontier Airlines Inc.
## 8 FL AirTran Airways Corporation
## 9 HA Hawaiian Airlines Inc.
## 10 MQ Envoy Air
## 11 OO SkyWest Airlines Inc.
## 12 UA United Air Lines Inc.
## 13 US US Airways Inc.
## 14 VX Virgin America
## 15 WN Southwest Airlines Co.
## 16 YV Mesa Airlines Inc.
airports
## # A tibble: 1,458 × 8
## faa name lat lon alt tz dst tzone
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 04G Lansdowne Airport 41.1 -80.6 1044 -5 A America/…
## 2 06A Moton Field Municipal Airport 32.5 -85.7 264 -6 A America/…
## 3 06C Schaumburg Regional 42.0 -88.1 801 -6 A America/…
## 4 06N Randall Airport 41.4 -74.4 523 -5 A America/…
## 5 09J Jekyll Island Airport 31.1 -81.4 11 -5 A America/…
## 6 0A9 Elizabethton Municipal Airport 36.4 -82.2 1593 -5 A America/…
## 7 0G6 Williams County Airport 41.5 -84.5 730 -5 A America/…
## 8 0G7 Finger Lakes Regional Airport 42.9 -76.8 492 -5 A America/…
## 9 0P2 Shoestring Aviation Airfield 39.8 -76.6 1000 -5 U America/…
## 10 0S9 Jefferson County Intl 48.1 -123. 108 -8 A America/…
## # ℹ 1,448 more rows
planes
## # A tibble: 3,322 × 9
## tailnum year type manufacturer model engines seats speed engine
## <chr> <int> <chr> <chr> <chr> <int> <int> <int> <chr>
## 1 N10156 2004 Fixed wing multi… EMBRAER EMB-… 2 55 NA Turbo…
## 2 N102UW 1998 Fixed wing multi… AIRBUS INDU… A320… 2 182 NA Turbo…
## 3 N103US 1999 Fixed wing multi… AIRBUS INDU… A320… 2 182 NA Turbo…
## 4 N104UW 1999 Fixed wing multi… AIRBUS INDU… A320… 2 182 NA Turbo…
## 5 N10575 2002 Fixed wing multi… EMBRAER EMB-… 2 55 NA Turbo…
## 6 N105UW 1999 Fixed wing multi… AIRBUS INDU… A320… 2 182 NA Turbo…
## 7 N107US 1999 Fixed wing multi… AIRBUS INDU… A320… 2 182 NA Turbo…
## 8 N108UW 1999 Fixed wing multi… AIRBUS INDU… A320… 2 182 NA Turbo…
## 9 N109UW 1999 Fixed wing multi… AIRBUS INDU… A320… 2 182 NA Turbo…
## 10 N110UW 1999 Fixed wing multi… AIRBUS INDU… A320… 2 182 NA Turbo…
## # ℹ 3,312 more rows
weather
## # A tibble: 26,115 × 15
## origin year month day hour temp dewp humid wind_dir wind_speed
## <chr> <int> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 EWR 2013 1 1 1 39.0 26.1 59.4 270 10.4
## 2 EWR 2013 1 1 2 39.0 27.0 61.6 250 8.06
## 3 EWR 2013 1 1 3 39.0 28.0 64.4 240 11.5
## 4 EWR 2013 1 1 4 39.9 28.0 62.2 250 12.7
## 5 EWR 2013 1 1 5 39.0 28.0 64.4 260 12.7
## 6 EWR 2013 1 1 6 37.9 28.0 67.2 240 11.5
## 7 EWR 2013 1 1 7 39.0 28.0 64.4 240 15.0
## 8 EWR 2013 1 1 8 39.9 28.0 62.2 250 10.4
## 9 EWR 2013 1 1 9 39.9 28.0 62.2 260 15.0
## 10 EWR 2013 1 1 10 41 28.0 59.6 260 13.8
## # ℹ 26,105 more rows
## # ℹ 5 more variables: wind_gust <dbl>, precip <dbl>, pressure <dbl>,
## # visib <dbl>, time_hour <dttm>
primary keyが本当に観測を一意に識別しているか検証する。primary
keyが無い場合はmutate()とrow_number()等で代替keyを作る。
planes %>% count(tailnum) %>% filter(n > 1) # primary keyのcount()で、各エントリについてnが1より大きくないか調べる
## # A tibble: 0 × 2
## # ℹ 2 variables: tailnum <chr>, n <int>
# weatherの観測を一意に識別するには、5つの変数が必要, countは自動的にgroup_by()することに注意
weather %>% count(year, month, day, hour, origin) %>% filter(n > 1)
## # A tibble: 3 × 6
## year month day hour origin n
## <int> <int> <int> <int> <chr> <int>
## 1 2013 11 3 1 EWR 2
## 2 2013 11 3 1 JFK 2
## 3 2013 11 3 1 LGA 2
更新ジョイン mutating joinは2つの表(X, Y)の変数を組み合わせる。
keyとマッチする観測を探し出し、変数を1つの表からもう1つの表にコピーする。mutate()と同様、ジョイン関数は変数を右側に追加する。
flights2 <- flights %>% select(year:day, hour, origin, dest, tailnum, carrier)
flights2
## # A tibble: 336,776 × 8
## year month day hour origin dest tailnum carrier
## <int> <int> <int> <dbl> <chr> <chr> <chr> <chr>
## 1 2013 1 1 5 EWR IAH N14228 UA
## 2 2013 1 1 5 LGA IAH N24211 UA
## 3 2013 1 1 5 JFK MIA N619AA AA
## 4 2013 1 1 5 JFK BQN N804JB B6
## 5 2013 1 1 6 LGA ATL N668DN DL
## 6 2013 1 1 5 EWR ORD N39463 UA
## 7 2013 1 1 6 EWR FLL N516JB B6
## 8 2013 1 1 6 LGA IAD N829AS EV
## 9 2013 1 1 6 JFK MCO N593JB B6
## 10 2013 1 1 6 LGA ORD N3ALAA AA
## # ℹ 336,766 more rows
# flightsに航空会社の正式名を追加するため、airlinesとflightsをleft_joinで組み合わせる
flights2 %>% select(-origin, -dest) %>% left_join(airlines, by = "carrier")
## # A tibble: 336,776 × 7
## year month day hour tailnum carrier name
## <int> <int> <int> <dbl> <chr> <chr> <chr>
## 1 2013 1 1 5 N14228 UA United Air Lines Inc.
## 2 2013 1 1 5 N24211 UA United Air Lines Inc.
## 3 2013 1 1 5 N619AA AA American Airlines Inc.
## 4 2013 1 1 5 N804JB B6 JetBlue Airways
## 5 2013 1 1 6 N668DN DL Delta Air Lines Inc.
## 6 2013 1 1 5 N39463 UA United Air Lines Inc.
## 7 2013 1 1 6 N516JB B6 JetBlue Airways
## 8 2013 1 1 6 N829AS EV ExpressJet Airlines Inc.
## 9 2013 1 1 6 N593JB B6 JetBlue Airways
## 10 2013 1 1 6 N3ALAA AA American Airlines Inc.
## # ℹ 336,766 more rows
内部ジョイン inner join, keyが共通した行が抽出される、 X AND Y。一般に観測が失われるため分析には不向き。
x <- tribble(
~key, ~val_x,
1, "x1",
2, "x2",
3, "x3"
)
y <- tribble(
~key, ~val_y,
1, "y1",
2, "y2",
4, "y3"
)
x %>%
inner_join(y, by = "key")
## # A tibble: 2 × 3
## key val_x val_y
## <dbl> <chr> <chr>
## 1 1 x1 y1
## 2 2 x2 y2
内部ジョインが両方の表にある観測を保持する一方、外部ジョイン outer joinは、少なくとも1つの表にある観測を保持する。以下の3種類がある:
left join:Xの全観測を保持する、マッチが無くても元の観測が保持されるので一番良く使われるright join:Yの全観測を保持するfull join:XとY全ての観測を保持する、X OR Y重複キー
# 1つの表に重複キーがある場合
x <- tribble(
~key, ~val_x,
1, "x1",
2, "x2",
2, "x3",
1, "x4"
)
y <- tribble(
~key, ~val_y,
1, "y1",
2, "y2"
)
left_join(x, y, by = "key")
## # A tibble: 4 × 3
## key val_x val_y
## <dbl> <chr> <chr>
## 1 1 x1 y1
## 2 2 x2 y2
## 3 2 x3 y2
## 4 1 x4 y1
# 両方の表に重複キーが含まれる場合、どちらの表でも観測を一意に識別できないので、普通はエラーになる。joinすると全ての組み合わせが得られる
x <- tribble(
~key, ~val_x,
1, "x1",
2, "x2",
2, "x3",
3, "x4"
)
y <- tribble(
~key, ~val_y,
1, "y1",
2, "y2",
2, "y3",
3, "y4"
)
left_join(x, y, by = "key")
## Warning in left_join(x, y, by = "key"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 2 of `x` matches multiple rows in `y`.
## ℹ Row 2 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.
## # A tibble: 6 × 3
## key val_x val_y
## <dbl> <chr> <chr>
## 1 1 x1 y1
## 2 2 x2 y2
## 3 2 x2 y3
## 4 2 x3 y2
## 5 2 x3 y3
## 6 3 x4 y4
Key列の指定
defaultではby =NULLで、両方の表に現れる全ての変数を使い、natural
joinする
flights2 %>%
left_join(weather)
## Joining with `by = join_by(year, month, day, hour, origin)`
## # A tibble: 336,776 × 18
## year month day hour origin dest tailnum carrier temp dewp humid
## <int> <int> <int> <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 2013 1 1 5 EWR IAH N14228 UA 39.0 28.0 64.4
## 2 2013 1 1 5 LGA IAH N24211 UA 39.9 25.0 54.8
## 3 2013 1 1 5 JFK MIA N619AA AA 39.0 27.0 61.6
## 4 2013 1 1 5 JFK BQN N804JB B6 39.0 27.0 61.6
## 5 2013 1 1 6 LGA ATL N668DN DL 39.9 25.0 54.8
## 6 2013 1 1 5 EWR ORD N39463 UA 39.0 28.0 64.4
## 7 2013 1 1 6 EWR FLL N516JB B6 37.9 28.0 67.2
## 8 2013 1 1 6 LGA IAD N829AS EV 39.9 25.0 54.8
## 9 2013 1 1 6 JFK MCO N593JB B6 37.9 27.0 64.3
## 10 2013 1 1 6 LGA ORD N3ALAA AA 39.9 25.0 54.8
## # ℹ 336,766 more rows
## # ℹ 7 more variables: wind_dir <dbl>, wind_speed <dbl>, wind_gust <dbl>,
## # precip <dbl>, pressure <dbl>, visib <dbl>, time_hour <dttm>
# flight2とplanesにはどちらもyearがあるが、それぞれ異なる意味なのでtailnumだけでjoinしたい
flights2 %>%
left_join(planes, by = "tailnum")
## # A tibble: 336,776 × 16
## year.x month day hour origin dest tailnum carrier year.y type
## <int> <int> <int> <dbl> <chr> <chr> <chr> <chr> <int> <chr>
## 1 2013 1 1 5 EWR IAH N14228 UA 1999 Fixed wing mult…
## 2 2013 1 1 5 LGA IAH N24211 UA 1998 Fixed wing mult…
## 3 2013 1 1 5 JFK MIA N619AA AA 1990 Fixed wing mult…
## 4 2013 1 1 5 JFK BQN N804JB B6 2012 Fixed wing mult…
## 5 2013 1 1 6 LGA ATL N668DN DL 1991 Fixed wing mult…
## 6 2013 1 1 5 EWR ORD N39463 UA 2012 Fixed wing mult…
## 7 2013 1 1 6 EWR FLL N516JB B6 2000 Fixed wing mult…
## 8 2013 1 1 6 LGA IAD N829AS EV 1998 Fixed wing mult…
## 9 2013 1 1 6 JFK MCO N593JB B6 2004 Fixed wing mult…
## 10 2013 1 1 6 LGA ORD N3ALAA AA NA <NA>
## # ℹ 336,766 more rows
## # ℹ 6 more variables: manufacturer <chr>, model <chr>, engines <int>,
## # seats <int>, speed <int>, engine <chr>
# この場合、出力ではyear.x、year.yと区別される
# c("a" = "b")とすると、表Xのaと表Yのbをマッチさせる
flights2 %>%
left_join(airports, c("dest" = "faa"))
## # A tibble: 336,776 × 15
## year month day hour origin dest tailnum carrier name lat lon alt
## <int> <int> <int> <dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 2013 1 1 5 EWR IAH N14228 UA Georg… 30.0 -95.3 97
## 2 2013 1 1 5 LGA IAH N24211 UA Georg… 30.0 -95.3 97
## 3 2013 1 1 5 JFK MIA N619AA AA Miami… 25.8 -80.3 8
## 4 2013 1 1 5 JFK BQN N804JB B6 <NA> NA NA NA
## 5 2013 1 1 6 LGA ATL N668DN DL Harts… 33.6 -84.4 1026
## 6 2013 1 1 5 EWR ORD N39463 UA Chica… 42.0 -87.9 668
## 7 2013 1 1 6 EWR FLL N516JB B6 Fort … 26.1 -80.2 9
## 8 2013 1 1 6 LGA IAD N829AS EV Washi… 38.9 -77.5 313
## 9 2013 1 1 6 JFK MCO N593JB B6 Orlan… 28.4 -81.3 96
## 10 2013 1 1 6 LGA ORD N3ALAA AA Chica… 42.0 -87.9 668
## # ℹ 336,766 more rows
## # ℹ 3 more variables: tz <dbl>, dst <chr>, tzone <chr>
flights2 %>%
left_join(airports, c("origin" = "faa"))
## # A tibble: 336,776 × 15
## year month day hour origin dest tailnum carrier name lat lon alt
## <int> <int> <int> <dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 2013 1 1 5 EWR IAH N14228 UA Newar… 40.7 -74.2 18
## 2 2013 1 1 5 LGA IAH N24211 UA La Gu… 40.8 -73.9 22
## 3 2013 1 1 5 JFK MIA N619AA AA John … 40.6 -73.8 13
## 4 2013 1 1 5 JFK BQN N804JB B6 John … 40.6 -73.8 13
## 5 2013 1 1 6 LGA ATL N668DN DL La Gu… 40.8 -73.9 22
## 6 2013 1 1 5 EWR ORD N39463 UA Newar… 40.7 -74.2 18
## 7 2013 1 1 6 EWR FLL N516JB B6 Newar… 40.7 -74.2 18
## 8 2013 1 1 6 LGA IAD N829AS EV La Gu… 40.8 -73.9 22
## 9 2013 1 1 6 JFK MCO N593JB B6 John … 40.6 -73.8 13
## 10 2013 1 1 6 LGA ORD N3ALAA AA La Gu… 40.8 -73.9 22
## # ℹ 336,766 more rows
## # ℹ 3 more variables: tz <dbl>, dst <chr>, tzone <chr>
Filtering join フィルタジョイン、変数ではなく観測に影響を及ぼす。
semi_join(x, y):セミジョインはyのとマッチするxの全ての観測を保持する。anti_join(x, y):アンチジョインはyのとマッチするxの全ての観測を取り除く。# 例) 人々が最もよく行く旅行先のトップ10を探す。
top_dest <- flights %>% count(dest, sort = T) %>% head(10)
top_dest
## # A tibble: 10 × 2
## dest n
## <chr> <int>
## 1 ORD 17283
## 2 ATL 17215
## 3 LAX 16174
## 4 BOS 15508
## 5 MCO 14082
## 6 CLT 14064
## 7 SFO 13331
## 8 FLL 12055
## 9 MIA 11728
## 10 DCA 9705
# これらへのフライトを探す
flights %>% filter(dest %in% top_dest$dest)
## # A tibble: 141,145 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 542 540 2 923 850
## 2 2013 1 1 554 600 -6 812 837
## 3 2013 1 1 554 558 -4 740 728
## 4 2013 1 1 555 600 -5 913 854
## 5 2013 1 1 557 600 -3 838 846
## 6 2013 1 1 558 600 -2 753 745
## 7 2013 1 1 558 600 -2 924 917
## 8 2013 1 1 558 600 -2 923 937
## 9 2013 1 1 559 559 0 702 706
## 10 2013 1 1 600 600 0 851 858
## # ℹ 141,135 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
# しかしこの方式を複数の変数について拡張するのは容易ではないので、代わりにsemi_joinを使用する
flights %>% semi_join(top_dest)
## Joining with `by = join_by(dest)`
## # A tibble: 141,145 × 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 542 540 2 923 850
## 2 2013 1 1 554 600 -6 812 837
## 3 2013 1 1 554 558 -4 740 728
## 4 2013 1 1 555 600 -5 913 854
## 5 2013 1 1 557 600 -3 838 846
## 6 2013 1 1 558 600 -2 753 745
## 7 2013 1 1 558 600 -2 924 917
## 8 2013 1 1 558 600 -2 923 937
## 9 2013 1 1 559 559 0 702 706
## 10 2013 1 1 600 600 0 851 858
## # ℹ 141,135 more rows
## # ℹ 11 more variables: arr_delay <dbl>, carrier <chr>, flight <int>,
## # tailnum <chr>, origin <chr>, dest <chr>, air_time <dbl>, distance <dbl>,
## # hour <dbl>, minute <dbl>, time_hour <dttm>
アンチジョインはジョインのミスマッチの診断に使うことができる。
# flightsとplanesを連結する時に、planesとマッチがとれないflightsが多いかどうか調べる
flights %>% anti_join(planes, by = "tailnum") %>% count(tailnum, sort = T)
## # A tibble: 722 × 2
## tailnum n
## <chr> <int>
## 1 <NA> 2512
## 2 N725MQ 575
## 3 N722MQ 513
## 4 N723MQ 507
## 5 N713MQ 483
## 6 N735MQ 396
## 7 N0EGMQ 371
## 8 N534MQ 364
## 9 N542MQ 363
## 10 N531MQ 349
## # ℹ 712 more rows
intersect(x,y):xとyに共通な観測だけを返す。
union(x,y):xとyにある観測を一意にして返す。
setdiff(x,y):xにはあるがyにはない観測を返す。
df1 <- tribble(
~x, ~y,
1, 1,
2, 1
)
df2 <- tribble(
~x, ~y,
1, 1,
1, 2
)
intersect(df1, df2)
## # A tibble: 1 × 2
## x y
## <dbl> <dbl>
## 1 1 1
union(df1, df2)
## # A tibble: 3 × 2
## x y
## <dbl> <dbl>
## 1 1 1
## 2 2 1
## 3 1 2
setdiff(df1, df2)
## # A tibble: 1 × 2
## x y
## <dbl> <dbl>
## 1 2 1
setdiff(df2, df1)
## # A tibble: 1 × 2
## x y
## <dbl> <dbl>
## 1 1 2
tibbleにはデータフレームやtibbleも格納できる。tidyr::nest()を使うと、データをtibble内に畳み込むことができる。畳み込まれたサブデータはデフォルトではリスト列”data”に格納される。名前を変更する場合は引数.keyで指定する。畳み込まれたデータ構造のことをネストデータあるいは入れ子データと呼ぶ。
# ダミーデータの生成, LETTERSとlettersは共にbuilt-in constantsでそれぞれ大文字と小文字のアルファベットが格納されている
set.seed(71)
N <- 15
dat <- tibble(tag1 = sample(LETTERS[1:3], N, replace = TRUE),
tag2 = sample(letters[1:5], N, replace = TRUE),
y = rnorm(N),
x = runif(N))
tibbleを要素に持つ、リスト形式で格納されている。
dat %>%
nest(.key = "data")
## # A tibble: 1 × 1
## data
## <list>
## 1 <tibble [15 × 4]>
事前に変数tag1でグループ化しておくと、tag1の水準ごとにデータを畳み込む。
dat %>%
group_by(tag1) %>%
nest()
## # A tibble: 3 × 2
## # Groups: tag1 [3]
## tag1 data
## <chr> <list>
## 1 C <tibble [6 × 3]>
## 2 B <tibble [5 × 3]>
## 3 A <tibble [4 × 3]>
group_nest()を使うとグループ化と畳み込みを同時に行う。
dat %>%
group_nest(tag1, tag2)
## # A tibble: 9 × 3
## tag1 tag2 data
## <chr> <chr> <list<tibble[,2]>>
## 1 A a [1 × 2]
## 2 A c [1 × 2]
## 3 A d [2 × 2]
## 4 B b [1 × 2]
## 5 B c [2 × 2]
## 6 B e [2 × 2]
## 7 C b [2 × 2]
## 8 C d [1 × 2]
## 9 C e [3 × 2]
畳み込まれたデータを広げるにはunnest()を使う。リスト列が複数ある場合、広げるリスト列を()内で指定できる。
dat %>%
group_nest(tag1, tag2) %>%
unnest()
## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(data)`.
## # A tibble: 15 × 4
## tag1 tag2 y x
## <chr> <chr> <dbl> <dbl>
## 1 A a 0.324 0.0138
## 2 A c -0.904 0.414
## 3 A d 0.900 0.427
## 4 A d 1.29 0.555
## 5 B b -1.62 0.473
## 6 B c -0.730 0.161
## 7 B c 1.70 0.889
## 8 B e 2.37 0.239
## 9 B e -0.673 0.648
## 10 C b 0.252 0.550
## 11 C b 0.324 0.397
## 12 C d -1.15 0.449
## 13 C e 0.0100 0.808
## 14 C e -0.887 0.0279
## 15 C e -1.79 0.110
このような入れ子データ構造は、グループ分けされたサブデータに対してmap関数で統計処理を適用するのに有用。
data列の各データのxの平均値を計算して新しい列に格納してみる。map()だとリストが返り値になるが、map_dbl()だと実数が返り値。
dat_nest <- dat %>%
group_nest(tag1)
dat_nest %>%
mutate(x_mean = map_dbl(data, ~ mean(.$x))) # .でリスト列中の各リストを表すためには~による記法が必要
## # A tibble: 3 × 3
## tag1 data x_mean
## <chr> <list<tibble[,3]>> <dbl>
## 1 A [4 × 3] 0.352
## 2 B [5 × 3] 0.482
## 3 C [6 × 3] 0.390
group_by()とsummarise()を使用した場合と比較して、統計要約量を別のデータフレームに格納する必要が無いのが利点。同様に統計検定も可能。
dat_nest %>%
mutate(ttest = map(data, ~ t.test(.$x, .$y)),
pval = map_dbl(ttest, ~.$p.value))
## # A tibble: 3 × 4
## tag1 data ttest pval
## <chr> <list<tibble[,3]>> <list> <dbl>
## 1 A [4 × 3] <htest> 0.926
## 2 B [5 × 3] <htest> 0.745
## 3 C [6 × 3] <htest> 0.0456
各サブデータのさらに一部を抽出して検定することも可能。以下ではtag2 == "b"でないデータについて検定している。
dat_nest %>%
mutate(x = map(data, ~ filter(., tag2 != "b") %>% .$x),
y = map(data, ~ filter(., tag2 != "b") %>% .$y),
ttest = map2(x, y, ~ t.test(.x, .y)),
pval = map_dbl(ttest, ~ .$p.value))
## # A tibble: 3 × 6
## tag1 data x y ttest pval
## <chr> <list<tibble[,3]>> <list> <list> <list> <dbl>
## 1 A [4 × 3] <dbl [4]> <dbl [4]> <htest> 0.926
## 2 B [5 × 3] <dbl [4]> <dbl [4]> <htest> 0.837
## 3 C [6 × 3] <dbl [4]> <dbl [4]> <htest> 0.0311
ggplot2はHadely Wickhamによって開発されたデータ可視化のためのパッケージ。「Grammer of Graphics」の思想に基づいている。グラフの構成要素のそれぞれを1つの層(レイヤー)と捉え、これらを重ねることでグラフを作成する。
geom_line()といったように、幾何オブジェクトはgeom_で始まる関数で表される。mapping = aes(マッピング)で指定する。mapping=は省略可能。まずデータの入ったデータフレームを呼び出し、ggplot()関数でデフォルトのマッピングを指定する。下の例ではxに変数coutryを、yに変数casesを指定。次にどのような幾何オブジェクトでプロットするかを指定する。下の例ではgeom_point()で散布図にする。またgeom_point()内にaes()を置き、変数yearを点の色であるcolorにマッピングする。このcolor = yearというマッピングはgeom_point()に限定されており、他の幾何オブジェクトには作用しない。もしggplot()内のaes()でcolor = yearとマッピングすれば、以下の全ての幾何オブジェクトに作用する。
library(tidyverse)
tb1<- table1 # tidyverseに収録されているデータ
plot1 <- tb1 %>%
ggplot(mapping = aes(x = country, y = cases)) +
geom_point(aes(color = year))
plot1
マッピングと設定の違い。colorの値をaes()の外でredに指定する。データの値とは無関係に全ての点が赤色になる。
plot2 <- tb1 %>%
ggplot(mapping = aes(x = country, y = cases)) +
geom_point(color = "red")
plot2
デフォルトでは変数yearは数値型であり連続的変数としてcolorにマッピングされた。factor型に変換してから離散的変数としてマッピングすることもできる。
plot3 <- tb1 %>%
mutate(year = as.factor(year)) %>%
ggplot(mapping = aes(x = country, y = cases)) +
geom_point(aes(color = year))
plot3
スケール変更の例。scale_color_manual()で色のスケール(色の対応のさせ方)を変更。またscale_y_continuous()でy軸の範囲を変更する。
plot3r <- tb1 %>%
mutate(year = as.factor(year)) %>%
ggplot(mapping = aes(x = country, y = cases)) +
geom_point(aes(color = year)) +
scale_colour_manual(values = c("orange", "forestgreen"))+
scale_y_continuous(limits = c(0, 250000))
plot3r
ggplot2において色はエステティック属性の1つであり、x座標やサイズ(size)と同じように扱われる。
幾何オブジェクトの色を設定するには、geom_**()関数のaes()の
外でcolorまたはfillの値を設定する。
変数を幾何オブジェクトの色にマッピングするには、aes()の中でcolorまたはfillの値に変数をマッピングする。
library(RColorBrewer)
display.brewer.all(type = "div")
display.brewer.all(type = "seq")
display.brewer.all(type = "qual")
display.brewer.all(colorblindFriendly = TRUE)
brewer.pal(4, "Set2")
## [1] "#66C2A5" "#FC8D62" "#8DA0CB" "#E78AC3"
dput(brewer.pal(4, "Set2")) # dput()はオブジェクトをテキスト形式でコンソールに出力する、ベクトルを生成するのに便利
## c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3")
display.brewer.pal(4, "Set2")
brewer.pal.info
## maxcolors category colorblind
## BrBG 11 div TRUE
## PiYG 11 div TRUE
## PRGn 11 div TRUE
## PuOr 11 div TRUE
## RdBu 11 div TRUE
## RdGy 11 div FALSE
## RdYlBu 11 div TRUE
## RdYlGn 11 div FALSE
## Spectral 11 div FALSE
## Accent 8 qual FALSE
## Dark2 8 qual TRUE
## Paired 12 qual TRUE
## Pastel1 9 qual FALSE
## Pastel2 8 qual FALSE
## Set1 9 qual FALSE
## Set2 8 qual TRUE
## Set3 12 qual FALSE
## Blues 9 seq TRUE
## BuGn 9 seq TRUE
## BuPu 9 seq TRUE
## GnBu 9 seq TRUE
## Greens 9 seq TRUE
## Greys 9 seq TRUE
## Oranges 9 seq TRUE
## OrRd 9 seq TRUE
## PuBu 9 seq TRUE
## PuBuGn 9 seq TRUE
## PuRd 9 seq TRUE
## Purples 9 seq TRUE
## RdPu 9 seq TRUE
## Reds 9 seq TRUE
## YlGn 9 seq TRUE
## YlGnBu 9 seq TRUE
## YlOrBr 9 seq TRUE
## YlOrRd 9 seq TRUE
ggplot2の標準のカラーパレット
dum_dat <- data.frame(key=c("A", "B", "C", "D", "E", "F", "G", "H"),
value = c(1, 1, 1, 1, 1, 1, 1, 1), stringsAsFactors = TRUE)
std_col_pal <- dum_dat %>%
ggplot(aes(x = key, y = value, fill=key))+
geom_col()+
scale_fill_manual(values = c("#F8766D", "#CD9600", "#7CAE00", "#00BE67", "#00BFC4", "#00A9FF", "#C77CFF", "#FF61CC"))+
xlab("")+
ylab("")
std_col_pal
無難な3色
acceptable_3col_pal <- dum_dat %>%
filter(key %in% c("A", "B", "C")) %>%
ggplot(aes(x = key, y = value, fill=key))+
geom_col()+
scale_fill_manual(values = c('#bdbdbd','#3182bd',"#e6550d"))+
xlab("")+
ylab("")
acceptable_3col_pal
無難な3色(やや明るめ)
acceptable_3col_pal <- dum_dat %>%
filter(key %in% c("A", "B", "C")) %>%
ggplot(aes(x = key, y = value, fill=key))+
geom_col()+
# scale_fill_manual(values = c("#c1afa8", "#7ec8ef","#F8766D"))+
# scale_fill_manual(values = c("#a9a9ac","#5da6f4", "#ff8078"))+
# scale_fill_manual(values = c("#c1afa8","#5da6f4", "#ff8078"))+
scale_fill_manual(values = c("#bdbdbd","#5da6f4", "#F8766D"))+
xlab("")+
ylab("")
acceptable_3col_pal
無難な5色
acceptable_5col_pal <- dum_dat %>%
filter(key %in% c("A", "B", "C", "D", "E")) %>%
ggplot(aes(x = key, y = value, fill=key))+
geom_col()+
scale_fill_manual(values = c('#bdbdbd',"#F8766D", "#00BFC4", "#619CFF", "#F564E3"))+
xlab("")+
ylab("")
acceptable_5col_pal
モデルの目的は、データセットの簡単な低次元の要約。仮説生成と仮説確認には別のデータセットを使用する。
1. モデルのfamilyを定義する
2. 適合モデル(a fitted
model)をデータに最も近いファミリーから探索し生成する、適合モデルは何らかの基準で「最良」のモデル
パッケージmodelrは基本Rのモデル化機能をパイプで自然に使えるようwrapしたもの。
library(tidyverse)
library(modelr)
options(na.action = na.warn)
sim1はmodelrに格納されている模擬データセットの1つ、線形パターンy = a_0 + a_1 * xを示す。
ggplot(sim1, aes(x,y)) +
geom_point()
sim1
## # A tibble: 30 × 2
## x y
## <int> <dbl>
## 1 1 4.20
## 2 1 7.51
## 3 1 2.13
## 4 2 8.99
## 5 2 10.2
## 6 2 11.3
## 7 3 7.36
## 8 3 10.5
## 9 3 10.5
## 10 4 12.4
## # ℹ 20 more rows
ランダムに線を重ねる。
models <- tibble(
# -20〜40の間で250個の実数をランダムに生成
a1 = runif(250, -20, 40),
a2 = runif(250, -5, 5)
)
ggplot(sim1, aes(x, y)) +
# 250組の(a1, a2)により250本の直線が引かれる
geom_abline(aes(intercept = a1, slope = a2), data = models, alpha = 1/4) +
geom_point()
距離計算
model1はモデルのファミリーを表す、モデルの引数(xの値)とデータを入力にとり、モデルの予測値を出力。
model1 <- function(a, data){
a[1] + data$x * a[2] ##data$xはデータフレームdataの列xを表す
}
model1(c(7, 1.5), sim1)
## [1] 8.5 8.5 8.5 10.0 10.0 10.0 11.5 11.5 11.5 13.0 13.0 13.0 14.5 14.5 14.5
## [16] 16.0 16.0 16.0 17.5 17.5 17.5 19.0 19.0 19.0 20.5 20.5 20.5 22.0 22.0 22.0
モデルによるy値(予測)と実際のデータのy値(応答)の差 = 距離を計算する。複数の個別の距離を全体的な距離にまとめる必要がある、これには「偏差の二乗平均平方根 root-mean squared deviation」(RMSD)を使う。実際と予測の差を計算して、二乗し、平均を取り、その平方根を求める。
measure_distance <- function(mod, data){
diff <- data$y - model1(mod, data)
sqrt(mean(diff ^ 2))
}
measure_distance(c(7, 1.5), sim1)
## [1] 2.665212
modelsの中の250組全ての(a1, a2)に対してRMSDを計算。
sim1_dist <- function(a1, a2){
measure_distance(c(a1, a2), sim1)
}
# map2は引数を2つ取ることができる、計算したRMSDを列distに格納
models <- models %>%
mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
models
## # A tibble: 250 × 3
## a1 a2 dist
## <dbl> <dbl> <dbl>
## 1 -10.6 -1.67 36.9
## 2 -15.4 -0.839 36.5
## 3 9.19 -2.13 21.8
## 4 -3.37 4.42 8.94
## 5 -13.6 -1.86 40.9
## 6 -14.6 2.27 17.8
## 7 38.9 -4.28 18.3
## 8 30.1 1.15 21.1
## 9 33.5 0.688 22.2
## 10 1.85 -2.52 30.6
## # ℹ 240 more rows
最良の10モデルをプロット、distが小さいほど良いモデル。
ggplot(sim1, aes(x, y)) +
geom_point(size = 2, colour = "grey30") +
geom_abline(
aes(intercept = a1, slope = a2, colour = -dist), ## color = -distで最小のdistが一番明るい色になる
data = filter(models, rank(dist) <= 10) ## distの小さいもの10位を選ぶ
)
モデルを観測と考えて、散布図を描き, -distで色付け、最良の10モデルを赤丸で囲む。
ggplot(models, aes(a1, a2)) +
geom_point(
data = filter(models, rank(dist) <= 10),
size =4, color = "red"
) +
geom_point(aes(color = -dist))
格子探索、より系統的に等間隔に格子点を生成する。どこに最良のモデルがあるか見当がついたので、大雑把に格子のパラメータを選ぶ。expand.grid関数はn個のベクトルに含まれる要素の全ての組み合わせをすばやく書き出す。
grid <- expand.grid(
a1 = seq(-5, 20, length = 25),
a2 = seq(1, 3, length = 25)
) %>%
mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
grid %>%
ggplot(aes(a1, a2)) +
geom_point(
data = filter(grid, rank(dist) <= 10),
size = 4, color = "red"
) +
geom_point(aes(color = -dist))
最良の10モデルを元のデータに重ねる。
ggplot(sim1, aes(x, y)) +
geom_point(size = 2, color = "grey30") +
geom_abline(
aes(intercept = a1, slope = a2, color = -dist),
data = filter(grid, rank(dist) <= 10)
)
ニュートン法による最良モデルの探索
モデルとデータセットとの距離を定義する関数とモデルのパラメータを変更して距離を最小にできるアルゴリズムがあれば、最良のモデルが求まる。
best <- optim(c(0,0), measure_distance, data = sim1)
best$par
## [1] 4.222248 2.051204
ggplot(sim1, aes(x,y)) +
geom_point(size =2, color = "grey30") +
geom_abline(intercept = best$par[1], slope = best$par[2])
線形モデルの適合専用に設計されたlm()はoptim()を使わないで、線形モデルの数学的構造を活用する。lm()はフォーミュラ(formula)という特別な形式でモデルファミリーを指定する。y ~ xはy = a_1 + a_2 * xと同等。
sim1_mod <- lm(y ~ x, data = sim1) ##
coef(sim1_mod)
## (Intercept) x
## 4.220822 2.051533
予測はモデルが補足したパターンを示唆。
data_grid()はデータのある領域を等間隔に覆う格子を生成。modeler::add_prediction()は引数にデータフレームとモデルを取り、モデルの予測をデータフレームの新たな列に追加する。
grid <- sim1 %>%
data_grid(x)
grid
## # A tibble: 10 × 1
## x
## <int>
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
## 10 10
#予測を追加、列predが追加される
grid <- grid %>%
add_predictions(sim1_mod)
grid
## # A tibble: 10 × 2
## x pred
## <int> <dbl>
## 1 1 6.27
## 2 2 8.32
## 3 3 10.4
## 4 4 12.4
## 5 5 14.5
## 6 6 16.5
## 7 7 18.6
## 8 8 20.6
## 9 9 22.7
## 10 10 24.7
ggplot(sim1, aes(x)) + ##予測をプロット
geom_point(aes(y=y)) +
geom_line(
aes(y = pred),
data= grid,
colour = "red",
linewidth = 1
)
残差とは観測値と予測値の距離。残差はモデルが見逃しているパターンを伝える。modeler::add_residuals()は引数に元のデータセットとモデルを取り、データフレームに残差を追加。残差の平均は常に0になる。
sim1 <- sim1 %>%
add_residuals(sim1_mod) ##残差 residが追加される
sim1
## # A tibble: 30 × 3
## x y resid
## <int> <dbl> <dbl>
## 1 1 4.20 -2.07
## 2 1 7.51 1.24
## 3 1 2.13 -4.15
## 4 2 8.99 0.665
## 5 2 10.2 1.92
## 6 2 11.3 2.97
## 7 3 7.36 -3.02
## 8 3 10.5 0.130
## 9 3 10.5 0.136
## 10 4 12.4 0.00763
## # ℹ 20 more rows
# 度数分布多角形を描いて残差の分布の広がりを理解する
ggplot(sim1, aes(resid)) +
geom_freqpoly(binwidth = 0.5)
元の予測の代わりに残差を使ってプロットすることも多い。残差がランダムなノイズに見えるなら、モデルはデータセットのパターンをしっかり補足したと言える。
ggplot(sim1, aes(x, resid)) +
geom_ref_line(h = 0) + ## 参照となる線を引く
geom_point()
model_matrix()はデータフレームとフォーミュラを引数にし、モデル方程式を表すtibbleを返す。出力の各列にはモデルの係数が関連づけられ、関数は常にy = a_1 * out1 + a_2 * out_2。
df <- tribble(
~y, ~x1, ~x2,
4, 2, 5,
5, 1, 6
)
# Rが切片をモデルに追加する時は、デフォルトで1ばかりの列を追加する
model_matrix(df, y ~ x1)
## # A tibble: 2 × 2
## `(Intercept)` x1
## <dbl> <dbl>
## 1 1 2
## 2 1 1
# 1ばかりの列を追加したくなければ、-1で明示的に取り消す
model_matrix(df, y ~ x1 - 1)
## # A tibble: 2 × 1
## x1
## <dbl>
## 1 2
## 2 1
model_matrix(df, y ~ x1 + x2)
## # A tibble: 2 × 3
## `(Intercept)` x1 x2
## <dbl> <dbl> <dbl>
## 1 1 2 5
## 2 1 1 6
model_matrix()はカテゴリ変数を(0,
1)の数値に変換できる。下記の例ではカテゴリ変数sexをsex_maleに変換。sex_maleはsexがmaleならば1、その他ならば0。
df <- tribble(
~ sex, ~ response,
"male", 1,
"female", 2,
"male", 1
)
model_matrix(df, response ~ sex)
## # A tibble: 3 × 2
## `(Intercept)` sexmale
## <dbl> <dbl>
## 1 1 1
## 2 1 0
## 3 1 1
sim2はカテゴリ変数x(a, b, c,
d)を含むサンプルデータ。カテゴリ変数xのモデルは、各カテゴリで平均値を予測する。平均が偏差の二乗平均平方根距離(RMSD)を最小化するから。
data(sim2)
ggplot(sim2) +
geom_point(aes(x,y))
mod2 <- lm(y ~ x, data = sim2) ## モデルを適合させて予測を生成
grid <- sim2 %>%
data_grid(x) %>%
add_predictions(mod2)
grid
## # A tibble: 4 × 2
## x pred
## <chr> <dbl>
## 1 a 1.15
## 2 b 8.12
## 3 c 6.13
## 4 d 1.91
ggplot(sim2, aes(x)) +
geom_point(aes(y = y)) +
geom_point(
data = grid,
aes(y = pred),
color = "red",
size =4
)
sim3にはカテゴリ予測子x2と連続予測子x1が含まれる。
data(sim3)
ggplot(sim3, aes(x1, y)) +
geom_point(aes(color = x2))
mod1:変数を+で追加するとモデルは効果を他の全てとは独立に推定。mod2:*で追加すると独立作用に加えて交互作用を考慮する。すなわちy ~ x1 * x2はy = a_0 + a_1 * x1 + a2_ * x2 + a_12 * x1 * x2、に翻訳される。
mod1 <- lm(y ~ x1 + x2, data = sim3) ## 変数を+で追加するとモデルは効果を他の全てとは独立に推定
mod2 <- lm(y ~ x1 * x2, data = sim3) ## *で追加すると独立作用に加えて交互作用を考慮する
2つの予測子があるのでdata_grid()に両方の変数を指定する必要がある、x1とx2の全ての一意な値を探し出し、全ての組み合わせを生成。両方のモデルから同時に予測を生成するには、spread_prediction()で各予測に対応する新たな列を作る。
grid <- sim3 %>%
data_grid(x1, x2) %>%
spread_predictions(mod1, mod2)
grid
## # A tibble: 40 × 4
## x1 x2 mod1 mod2
## <int> <fct> <dbl> <dbl>
## 1 1 a 1.67 1.21
## 2 1 b 4.56 7.52
## 3 1 c 6.48 5.71
## 4 1 d 4.03 2.32
## 5 2 a 1.48 1.12
## 6 2 b 4.37 6.66
## 7 2 c 6.28 5.68
## 8 2 d 3.84 2.50
## 9 3 a 1.28 1.02
## 10 3 b 4.17 5.81
## # ℹ 30 more rows
各予測を1列predにまとめて出力するには、行に追加するgather_prediction()を使用。
grid <- sim3 %>%
data_grid(x1, x2) %>%
gather_predictions(mod1, mod2)
grid
## # A tibble: 80 × 4
## model x1 x2 pred
## <chr> <int> <fct> <dbl>
## 1 mod1 1 a 1.67
## 2 mod1 1 b 4.56
## 3 mod1 1 c 6.48
## 4 mod1 1 d 4.03
## 5 mod1 2 a 1.48
## 6 mod1 2 b 4.37
## 7 mod1 2 c 6.28
## 8 mod1 2 d 3.84
## 9 mod1 3 a 1.28
## 10 mod1 3 b 4.17
## # ℹ 70 more rows
両方のモデルの結果をfacet_wrap()で1つのプロットにまとめる。2つの変数の独立作用のみを考慮したmod1ではどの線の傾きも等しい。交互作用を考慮したmod2では各線で傾きも切片も異なる。
ggplot(sim3, aes(x1, y, color = x2)) +
geom_point() +
geom_line(data = grid, aes(y = pred)) +
facet_wrap(~ model)
残差を調べて、どちらのモデルが適切かを検証。数学的手法ではなく、モデルが対象としているパターンを補足したかどうかの定性的評価に着目。mod2の残差では自明なパターンがほとんど無いのに対し、mod1の残差では、bである種のパターンが見逃されている。
sim3 <- sim3 %>%
gather_residuals(mod1, mod2)
ggplot(sim3, aes(x1, resid, color =x2)) +
geom_point() +
facet_grid(model~x2)
sim4は2つの連続変数x1, x2を含むサンプルデータ。交互作用を考慮しない場合とする場合でモデルを構築。
data(sim4) ## 2つの連続変数x1, x2を含む
mod1 <- lm(y ~ x1 + x2, data = sim4)
mod2 <- lm(y ~ x1 * x2, data = sim4)
grid <- sim4 %>%
data_grid(
x1 = seq_range(x1, 5), ## x1の一意の値すべてを使う代わりに最小値と最大値の間の5つの格子値を等間隔に使用
x2 = seq_range(x2, 5)
) %>%
gather_predictions(mod1, mod2)
grid
## # A tibble: 50 × 4
## model x1 x2 pred
## <chr> <dbl> <dbl> <dbl>
## 1 mod1 -1 -1 0.996
## 2 mod1 -1 -0.5 -0.395
## 3 mod1 -1 0 -1.79
## 4 mod1 -1 0.5 -3.18
## 5 mod1 -1 1 -4.57
## 6 mod1 -0.5 -1 1.91
## 7 mod1 -0.5 -0.5 0.516
## 8 mod1 -0.5 0 -0.875
## 9 mod1 -0.5 0.5 -2.27
## 10 mod1 -0.5 1 -3.66
## # ℹ 40 more rows
モデルを可視化、2つの連続予測子があるのでモデルを3次元表面で考える。
ggplot(grid, aes(x1, x2)) +
geom_tile(aes(fill = pred)) +
facet_wrap( ~ model)
表面を上から見るのではなく、複数のスライスを表示して横から見る。
ggplot(grid, aes(x1, pred, colour = x2, group = x2)) +
geom_line() +
facet_wrap(~ model)
ggplot(grid, aes(x2, pred, colour = x1, group = x1)) +
geom_line() +
facet_wrap(~ model)
log(y) ~ sqrt(x1) + sqrt(x2)はlog(y) = a_1 + a_2 * sqrt(x1) + a_3 * sqrt(x2)、に変換される。
変換が+,*,^,-を含むならI()でラップして、Rがモデル指定の一部と間違えないようにする。
例:
y ~ x + I(x^2)はy = a_1 + a_2*x + a_3 * x^2、と翻訳される。
I()を忘れて、y ~ x^2 + x、と指定するとRはy ~ x*x + xを計算する、x*xはxと自分自身の交互作用を意味するのでxと同じになり,
y = a_1 + a_2*xを指定することになる。
モデルが何をしているのかは、model_matrix()で方程式lm()が適合しているのかは正確に何であるかを確認できる。
df <- tribble(
~y, ~x,
1, 1,
2, 2,
3, 3
)
model_matrix(df, y ~ x^2 + x)
## # A tibble: 3 × 2
## `(Intercept)` x
## <dbl> <dbl>
## 1 1 1
## 2 1 2
## 3 1 3
model_matrix(df, y ~ I(x^2) + x)
## # A tibble: 3 × 3
## `(Intercept)` `I(x^2)` x
## <dbl> <dbl> <dbl>
## 1 1 1 1
## 2 1 4 2
## 3 1 9 3
poly()で多項式に適合させることができる。ただしデータの範囲外では、多項式が急激に正または負の無限大になるため、自然スプラインspline::ns()を使う方が安全。
model_matrix(df, y ~ poly(x, degree =2))
## # A tibble: 3 × 3
## `(Intercept)` `poly(x, degree = 2)1` `poly(x, degree = 2)2`
## <dbl> <dbl> <dbl>
## 1 1 -7.07e- 1 0.408
## 2 1 -9.07e-17 -0.816
## 3 1 7.07e- 1 0.408
library(splines)
model_matrix(df, y ~ ns(x, df = 2)) # dfはdegree of freedom
## # A tibble: 3 × 3
## `(Intercept)` `ns(x, df = 2)1` `ns(x, df = 2)2`
## <dbl> <dbl> <dbl>
## 1 1 0 0
## 2 1 0.566 -0.211
## 3 1 0.344 0.771
非線形関数の近似。
sim5 <- tibble(
x = seq(0, 3.5*pi, length = 50),
y = 4 * sin(x) + rnorm(length(x))
)
ggplot(sim5, aes(x,y)) + geom_point()
# 以下の5つのモデルを適合させる
mod1 <- lm(y ~ ns(x, 1), data = sim5)
mod2 <- lm(y ~ ns(x, 2), data = sim5)
mod3 <- lm(y ~ ns(x, 3), data = sim5)
mod4 <- lm(y ~ ns(x, 4), data = sim5)
mod5 <- lm(y ~ ns(x, 5), data = sim5)
grid <- sim5 %>%
data_grid(x = seq_range(x, n = 50, expand = 0.1)) %>% ## expand = 0.1は範囲を10%拡張する
gather_predictions(mod1, mod2, mod3, mod4, mod5, .pred = "y")
ggplot(sim5, aes(x, y)) +
geom_point() +
geom_line(data = grid, colour = "red") +
facet_wrap(~ model)
欠損値は変数間の関係に何の情報ももたらさないので、モデル関数は欠損値を含む行を削除する。
線形モデルでは残差が正規分布であることも仮定する。
一般化線形モデルglm()では目的変数の残差が正規分布に従わないデータに適用する、目的変数を変換して解析。
現実のデータを統計解析する際は、まずデータが正規分布に従っているか否か、すなわち正規性があるかどうかを確認する。
まずデータをプロットして分布の状態を目視で確認するのと同時に、正規性の検定を行う。
データの正規性の検定にはシャピロ・ウィルク検定がよく用いられる。Rではデフォルトでインストールされている関数shapiro.test()で検定することができる。この検定の帰無仮説は「データxが正規分布している」であり、計算されたp値が0.05以上ならば帰無仮説を棄却できず、従ってデータxは正規分布していると結論される。
データxが正規分布する例:
x <-rnorm(50, mean = 5, sd = 3)
shapiro.test(x)
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.9824, p-value = 0.6568
データxが正規分布しない例:
x <- runif(50, min = 0, max = 100)
shapiro.test(x)
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.94977, p-value = 0.03337
shapiro.test()に渡すデータにはNaNが含まれていても計算できるが、NaN以外の要素の個数が3〜5000という制限がある。そのため5000を超える場合はsample()による無作為抽出で対応する。
x <- rnorm(6000, mean = 10, sd = 3)
sample_5000 <- sample(x, 5000, replace = FALSE)
shapiro.test(sample_5000)
##
## Shapiro-Wilk normality test
##
## data: sample_5000
## W = 0.99962, p-value = 0.4695
shapiro.test()の入力はベクトルでありデータフレーム形式には対応していないため、データフレーム中のデータを渡すには自作関数を作る必要がある。実例はこちらのページ。
2群のデータを比較する場合、データに正規性がある場合はパラメトリック検定、正規性が無い場合はノンパラメトリック検定で統計検定する。
パラメトリック検定
ノンパラメトリック検定
3群以上のデータを比較する場合、データの正規性の検証後にさらに一元配置分散分析により各群の平均値に差があるか検証する必要がある。有意な差が検出された場合は、事後検定で比較する。
正規性がある場合
正規性が無い場合
各種検定の計算方法については、下記の実例集を参考にする。
下記のリンクを参照: